AN  INQUIRY:  EFFECTIVENESS  OF  THE  COMPLEX  EMPIRICAL  MODE 
DECOMPOSITION  METHOD,  THE  HILBERT-HUANG  TRANSFORM,  AND  THE 
FAS T-F OURIER  TRANSFORM  FOR  ANALYSIS  OF  DYNAMIC  OBJECTS 


THESIS 


Kristen  L.  Wallis,  Second  Lieutenant,  USAF 


AFIT/GE/ENG/12-42 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 
DISTRIBUTION  STATEMENT  A. 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the  United 
States  Government.  This  material  is  declared  a  work  of  the  U.S.  Government  and  is  not 
subject  to  copyright  protection  in  the  United  States. 


AFIT/GE/EN  G/ 1 2-42 


AN  INQUIRY:  EFFECTIVENESS  OF  THE  COMPLEX  EMPIRICAL  MODE 
DECOMPOSITION  METHOD,  THE  HILBERT-HUANG  TRANSFORM,  AND  THE 
FAS T-F OURIER  TRANSFORM  FOR  ANALYSIS  OF  DYNAMIC  OBJECTS 

THESIS 


Presented  to  the  Faculty 

Department  of  Electrical  and  Computer  Engineering 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Electrical  Engineering 


Kristen  L.  Wallis,  BSEE 
Second  Lieutenant,  USAF 
March  2012 

DISTRIBUTION  STATEMENT  A. 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AFIT/GE/EN  G/ 1 2-42 


AN  INQUIRY:  EFFECTIVENESS  OF  THE  COMPLEX  EMPIRICAL  MODE 
DECOMPOSITION  METHOD,  THE  HILBERT-HUANG  TRANSFORM,  AND  THE 
FAS T-F OURIER  TRANSFORM  FOR  ANALYSIS  OF  DYNAMIC  OBJECTS 


Kristen  L.  Wallis,  BSEE 
Second  Lieutenant,  USAF 


Approved: 


_ //signed// _ 

Andrew  J.  Terzuoli,  PhD  (Chairman) 


22  Feb  2012 
Date 


_ //signed// _ 

Lt  Col  Geoffrey  A.  Akers,  PhD  (Member) 


22  Feb  2012 
Date 


_ //signed// _ 

Peter  J.  Collins,  PhD  (Member) 


22  Feb  2012 
Date 


22  Feb  2012 


_ //signed// _ 

Mark  E.  Oxley,  PhD  (Member) 


Date 


AFIT/GE/EN  G/ 1 2-42 


Abstract 

A  review  of  current  signal  analysis  tools  show  that  new  techniques  are  required 
for  an  enhanced  fidelity  or  data  integrity.  Recently,  the  Hilbert-Huang  transform  (HHT) 
and  its  inherent  property,  the  Empirical  Mode  Decomposition  (EMD)  technique,  have 
been  formerly  investigated.  The  technique  of  Complex  EMD  (CEMD)  was  also 
explored.  The  scope  of  this  work  was  to  assess  the  CEMD  technique  as  an  innovative 
analysis  tool.  Subsequent  to  this,  comparisons  between  applications  of  the  Hilbert 
transform  (HT)  and  the  Fast-Fourier  transform  (FFT)  were  analyzed.  MATLAB®  was 
implemented  to  model  signal  decomposition  and  the  execution  of  mathematical 
transforms  for  generating  results.  The  CEMD  technique  successfully  decomposed  the 
data  into  its  oscillatory  modes.  After  comparative  graphical  analysis  of  the  HT  and  FFT, 
application  of  the  HT  provided  marginal  enhancements  of  the  data  modeled  previously  by 
the  FFT.  Altogether,  the  HHT  could  not  be  determined  as  a  helpful  analysis  tool. 
Nevertheless,  the  CEMD  technique,  an  inherent  component  of  the  HHT,  exhibited  a 
possible  improvement  as  an  analysis  tool  for  signal  processing  data.  Further  evaluation 
of  the  CEMD  technique  and  the  HHT  is  needed  for  ultimate  determination  of  their 
usefulness  as  an  analysis  tool. 
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AN  INQUIRY:  EFFECTIVENESS  OF  THE  COMPLEX  EMPIRICAL  MODE 
DECOMPOSITION  METHOD,  THE  HILBERT-HUANG  TRANSFORM,  AND  THE 
FAS T-F OURIER  TRANSFORM  FOR  ANALYSIS  OF  DYNAMIC  OBJECTS 


I.  Problem  Statement 


1.1.  Background 

The  effectiveness  of  the  Hilbert-Huang  transform  (HHT),  and  more  specifically 
the  complex  extension  of  the  Empirical  Mode  Decomposition  (EMD)  technique  within 
the  HHT,  as  analysis  tools  for  complex  radar-cross  section  (RCS)  data  collected  from 
dynamic  objects  were  being  investigated.  In  addition  to  assessing  the  use  of  the  HHT  for 
complex  RCS  data,  both  the  Fast-Fourier  transform  (FFT)  and  Hilbert  transform  were 
applied  to  the  data  that  has  been  analyzed  using  the  Complex  Empirical  Mode 
Decomposition  (CEMD)  algorithm.  Once  results  have  been  generated  for  both 
transforms  based  on  the  decomposed  complex  data,  the  outputs  from  the  two  transforms 
were  compared  to  determine  whether  the  HHT  can  provide  information  or  better  fidelity 
not  previously  provided  by  the  FFT  of  the  decomposed  complex  data. 

In  addition  to  assessing  the  EMD  technique  of  the  HHT,  the  EMD  technique  was 
extended  into  the  realm  of  complex-valued  signals  and  complex-valued  data  to  provide  a 
more  accurate  output  of  the  data  being  processed  by  the  FFT  and  the  Hilbert  transform  of 
the  decomposed  complex  data.  While  the  original  EMD  method  analyzes  only  the  real¬ 
valued  or  magnitude  portion  of  the  data,  the  use  of  an  algorithm  that  allows  the  analysis 
of  complex  data  was  implemented  as  part  of  assessing  and  examining  the  usefulness  of 
the  HHT. 


1 


1.1.1.  Benefits  of  Investigating  the  EMD  Technique 

The  benefits  of  the  findings  were  weighed  prior  to  beginning  the  research, 
especially  when  determining  their  validity  and  usefulness  as  a  possible  tool.  This 
particular  problem  does  not  require  the  use  of  expensive  equipment  or  many  undeveloped 
concepts,  but  there  is  a  cost  concerning  manpower  and  time  devoted  to  researching 
existing  topics  and  determining  how  the  new  topic  will  be  approached. 

For  the  task  of  assessing  the  CEMD  technique  and  HHT  as  analysis  tools,  as  well 
as  determining  the  usefulness  compared  with  the  FFT,  the  benefits  will  outweigh  the  time 
spent  in  the  research  and  development  towards  the  problem.  Due  to  the  previous  research 
done  in  the  field  of  knowledge  and  development  of  the  EMD  technique,  there  is  a  good 
foundation  for  the  development  of  the  CEMD  coding  that  will  be  necessary  for 
assessment  of  the  CEMD  method  and  the  use  of  the  HHT.  The  largest  benefit  of 
assessing  the  use  of  the  CEMD  method  and  the  HHT  presently  is  determining  its 
potential  for  future. 

1.1.2.  History  of  the  Hilbert- Huang  Transform 

Historically,  analysis  of  data  collected  from  dynamic  objects  has  been  done  using 
the  Fourier  spectral  analysis,  consisting  of  the  Discrete  Fourier  transform  (DFT)  and  the 
Fourier  transform.  In  the  mid-1960s,  after  the  development  of  the  FFT,  analysis  in 
computationally  intense  fields  was  revolutionized.  Even  though  the  Fourier  transform 
had  been  used  in  data  analysis  for  years  prior  to  the  introduction  of  the  FFT, 
implementation  of  the  FFT  algorithm  allowed  those  working  with  the  Fourier  spectral 
analysis  for  data  processing  experienced  an  expedited  form  of  signal  processing.  For  the 
past  50  years,  the  FFT  has  been  the  primary  analysis  tool  for  signal  processing.  Even 
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though  the  FFT  has  provided  an  acceptable  analysis  procedure  for  linear  and  stationary 
(or  periodic)  data  sets,  Fourier  analysis  has  some  constraints  and  drawbacks  inherent  to 
the  method. 

The  primary  drawback  inherent  to  Fourier  analysis,  and  by  extension  the  FFT 
algorithm,  is  a  limitation  founded  upon  the  analysis  of  linear  and  stationary  data.  The 
property  of  analysis  of  linear  and  stationary  data  is  due  to  the  FFT  employing  a  known 
and  predefined  set  of  basis  functions  that  maps  the  inputted  data  from  the  time  domain  to 
the  frequency  domain.  Consequently,  analysis  of  nonlinear  (such  as  shockwave  data  or 
turbulence  data)  and  non- stationary  data  is  nearly  impossible  with  the  FFT. 

As  a  result  of  the  constraint  on  the  FFT  method,  other  mathematical  tools  are 
receiving  consideration  from  analysts  in  the  intelligence  community.  In  recent  years,  an 
innovative  technique  for  analyzing  dynamic  objects  has  been  investigated  in  other  fields 
of  study.  The  specific  mathematical  tool  that  was  analyzed,  scrutinized,  and  assessed 
throughout  this  thesis  is  the  HHT,  a  more  recent  analysis  tool  that  contains  the  EMD 
technique  and  the  Hilbert  transform. 

Dr.  Huang  [1]  claimed  that  the  HHT  is  more  suitable  to  analyze  nonlinear  and 
non-stationary  objects,  resulting  in  an  interest  of  investigating  the  capabilities  and 
possibilities  of  using  the  HHT  as  an  analysis  tool.  The  HHT  performs  the  decomposition 
of  the  signal  into  its  oscillatory  modes  using  the  EMD  technique  and  applies  the  Hilbert 
transform  to  these  modes,  which  results  in  the  Hilbert  spectrum  representation  of  the 
signal. 
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1.1.3.  Explanation  of  the  Empirical  Mode  Decomposition  Technique 
The  EMD  algorithm  is  an  adaptive  technique,  in  which  a  given  signal  is 
decomposed  into  a  set  of  oscillating  basis  functions  through  use  of  the  sifting  process. 
These  oscillatory  components  of  the  decomposed  signal  are  called  the  Intrinsic  Mode 
Functions  (IMFs)  and  are  representations  of  the  oscillating  nature  embedded  in  the  data 
[1].  The  Fourier  transform  has  a  predefined  set  of  basis  functions,  yet  the  IMFs  that  are 
decomposed  from  the  data  are  considered  the  basis  functions  for  the  given  signal. 

More  precisely,  the  EMD  performs  the  following  decomposition  through  various 

steps: 


x{k)  =  '£jNj=ci{k)  +  r(k)  (1) 

where  ct(k),i  =  l,...,N  denote  the  IMFs  and  r(k )  denotes  the  residual.  The  IMF  is 

characterized  by  two  properties:  (1)  the  upper  and  lower  envelopes  are  symmetric;  and 
(2)  the  number  of  zero-crossings  and  the  number  of  extrema  are  exactly  equal  [2]. 

To  extract  the  IMFs  from  a  given  signal,  the  following  sifting  algorithm  is 
employed,  as  detailed  in  Table  1  [1], 

Table  1:  The  EMD  Algorithm 

1.  Find  the  locations  of  all  the  extrema  of  the  given  signal,  x  (k) 

2.  Interpolate  (using  the  cubic  spline  interpolation)  between  all  the  minima  and 
respective  maxima  to  obtain  the  signal  envelope  passing  through  the  minima 
emm(k)  and  respective  maxima  enrM (k ) 

3.  Compute  the  local  mean  m(k)  =  {enin{k)  +  enMX(k))l 2 

4.  Subtract  the  local  mean  from  the  signal  to  obtain  the  “oscillating”  signal 
s(k )  =  x\k)  —m(k) 

5.  If  the  resulting  signal  obeys  the  stopping  criterion,  d(k )  =  s(k)  becomes  an  IMF; 
otherwise,  set  x'(k)  =  s(k )  and  repeat  the  process  from  Step  1 
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The  stopping  criterion  for  the  final  step  is  the  normalized  squared  difference  between  two 


successive  iterates,  sn(k )  and  (k) ,  or  as  follows: 


jMK_,(*)-j,(*)ir  iSD 


i= 0 


sUk) 


(2) 


where  N  represents  the  total  number  of  samples  in  the  original  signal  x(k) ,  and  the 
standard  deviation  (SD)  is  set  within  the  range  of  (0.2-0. 3)  [1]. 

Once  the  sifting  algorithm  is  implemented,  the  Hilbert  transform  is  applied  to 
each  individual  IMF.  The  resulting  equation  upon  application  of  the  Hilbert  transform,  is 
given  by,  representing  the  Hilbert  spectrum  is  as  follows: 


X(t)  =  YJNl=a,(t)e 


(3) 


where  aft)  is  the  time-dependent  amplitude  and  Oft)  is  the  phase  function.  The 
instantaneous  radial  frequency  can  be  defined  by 


toft) 


dOft) 

clt 


(4) 


and  can  be  plotted  against  the  amplitude,  in  which  the  resulting  plot  is  the  Time- 
Frequency  Amplitude,  representing  the  Hilbert  spectrum.  The  combination  of  the 
instantaneous  frequency  concept,  detailed  in  (4),  and  the  EMD  technique,  detailed  in  (1), 
makes  the  HHT  a  powerful  tool  in  both  signal  decomposition  and  signal  analysis  [2]. 

1.1.4.  Primary  Focus  of  the  Hilbert-Huang  Transform 

The  HHT  consists  of  the  EMD  technique  and  Hilbert  transform,  in  which  the 
primary  focus  of  assessing  the  use  of  the  HHT  to  analyze  RCS  data  is  the  EMD  technique 
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portion.  The  Hilbert  spectrum  results  from  the  Hilbert  transform  application  and  is  used 
in  comparing  the  generated  graphical  results.  The  EMD  technique  has  been  investigated 
extensively  by  numerous  researchers  [3] -[8];  therefore,  for  purposes  of  data  analysis  in 
this  evaluation,  the  EMD  extension  into  the  complex  domain  will  be  explored  as  applied 
to  complex  RCS  data.  The  CEMD  method  has  become  more  relevant  in  fields  where  the 
data  collected  contains  a  phase  component  (such  as  signal  processing). 

1.1.5.  Complex  Extensions  of  the  Empirical  Mode  Decomposition  Technique 

Several  complex  extensions  of  the  EMD  technique  have  been  developed  recently. 
Such  complex  methods  include  CEMD  [9],  RICEMD  [10],  and  BEMD  [11].  Throughout 
literature,  the  term  CEMD  is  ambiguous,  referring  to  the  method  described  in  [9]  or  in 
reference  to  the  process  of  extending  the  EMD  technique  to  the  complex  domain. 

The  first  extension  introduced  for  extending  the  EMD  technique  in  the  complex 
domain  was  termed  CEMD  by  Tanaka  and  Mandic  [9].  The  authors  used  the  relationship 
that  exists  between  the  positive  and  negative  frequency  components  of  the  complex 
signal.  Rather  than  viewing  the  result  as  one  signal  consisting  of  these  two  frequency 
components,  Tanaka  and  Mandic  view  the  IMFs  extracted  as  independent  of  each  other. 
Then,  by  applying  the  EMD  technique  to  the  negative  and  positive  frequency 
components,  the  two  sets  of  IMFs  were  created,  which  corresponded  to  the  negative  and 
positive  frequency  components.  This  method  works  well  for  the  low-dimensional  (such  a 
single  dimensional)  case,  but  when  applied  to  a  higher-dimensional  case,  the  analysis 
degrades  due  to  the  rigorous  mathematical  nature  of  the  algorithm  and  the  non-intuitive 
extension  from  the  original  EMD  algorithm.  [9] 
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A  second  complex  EMD  extension  proposed  by  Umair  Bin  Altaf  et  al.  is 
RICEMD  [10].  This  method  differs  from  the  CEMD  method  developed  by  Tanaka  and 
Mandic  [9],  where  the  RICEMD  method  operates  completely  in  the  complex  domain.  As 
a  result,  the  steps  to  accomplish  the  EMD  method  are  performed  in  the  complex  domain 
rather  than  the  real  domain.  The  only  difference  between  the  original  EMD  method  and 
this  proposed  complex  extension  exists  in  determining  the  extrema  and  the  envelope  of 
the  signal.  Other  than  these  two  steps  of  the  EMD  algorithm  being  changed,  the 
RICEMD  and  the  original  EMD  methods  are  accomplished  in  the  same  manner. 

One  final  method  proposed  for  the  realization  of  the  complex  EMD  analysis  is 
called  BEMD  [11]  developed  by  Rilling  et  al.  This  method  is  based  on  the  idea  of 
“bivariate  signal  =  fast  rotations  +  slow  rotations.”  In  order  to  discriminate  between  the 
“fast”  and  “slow”  rotations,  the  idea  is  to  define  the  “slow”  component  as  the  mean  of  the 
defined  “envelope.”  The  envelope  for  the  bivariate  signal  is  now  three  dimensional 
rather  than  two  dimensional,  as  in  the  univariate  signal  case.  After  the  data  points  for 
analysis  of  the  signal  are  selected,  the  issue  of  defining  the  mean  arises.  The  preferred 
definition  for  the  mean  is  the  intersection  of  two  straight  lines,  where  the  lines 
intersecting  correspond  to  the  two  horizontal  tangents  and  the  two  vertical  tangents  [11]. 

The  desired  goal  for  defining  the  mean  was  the  same  as  the  original  EMD 
method:  a  smooth  curve  with  as  few  oscillations  as  possible.  The  interpolation  scheme 
employed  by  the  BEMD  technique  is  a  cubic  spline.  The  cubic  spline  was  employed  due 
to  its  minimum  curvature  property  and  fits  the  purpose  of  the  EMD  algorithm  best. 

The  BEMD  method  was  employed  to  extract  the  complex  IMFs  because  the  idea 
of  projecting  the  signal  into  specific  directions  to  extract  the  extrema,  as  well  as 
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connecting  the  points  to  create  the  desired  envelope  is  intuitive  and  similar  to  the  original 
EMD  method.  The  bivariate  time  series  in  the  algorithms  is  treated  as  a  complex-valued 
time  series.  Given  an  angle  direction,  the  bivariate  extensions  are  defined  by  the  EMD 
algorithm,  only  with  new  sifting  elementary  operators  defined  by  Sm  and  S112  .which 
correspond  to  the  two  algorithms  [11]. 

For  the  purposes  of  this  thesis,  the  BEMD  method  proposed  by  Rilling  et  al.  was 
employed  to  analyze  the  complex  RCS  data  provided  by  the  sponsor,  as  well  as  the 
simulated  data  collected  at  the  RCS  range. 

1.2.  Statement  of  Problem 

The  problem  addressed  in  this  research  effort  consisted  of  two  parts.  The  first 
part  was  assessing  the  use  of  the  CEMD  technique  of  the  HHT  as  an  analysis  tool  on 
complex  RCS  data  created  by  dynamic  objects.  Next,  the  decomposed  signal  after 
implementation  of  the  CEMD  method  was  analyzed  after  the  applications  of  the  FFT  and 
Hilbert  transform  to  determine  whether  the  HHT  provided  enhanced  analysis  of  the  signal 
as  compared  with  analysis  provided  previously  from  the  FFT. 

1.3.  Justification  for  Research 

Traditional  data  analysis  tools,  such  as  the  Fourier  transform  and  wavelet 
analysis,  depend  on  the  mapping  of  a  pre-existing  and  known  signal  being  transformed 
from  the  time  domain  to  the  frequency  domain  or  from  the  scale  domain  to  the  dilation 
domain,  respectively.  Such  tools  are  limited  to  the  analysis  of  linear  and  stationary 
signals,  leading  to  problems  when  a  need  to  analyze  nonlinear  and  non- stationary  data 
arises.  Most  data  collected  fall  within  this  category  of  nonlinear  and  non-stationary  data, 
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for  the  sake  of  simplicity,  there  is  a  standing  assumption  that  the  data  are  linear  and 
stationary.  Recently,  there  has  been  interest  in  analyzing  the  data  collected  from  dynamic 
objects  using  the  CEMD  technique  and  the  HHT.  Because  the  data  are  nonlinear  and 
non-stationary,  applying  this  recent  adaptive  data-driven  analysis  tool  may  provide  novel 
insights  about  the  data.  The  assessment  of  implementing  the  CEMD  method  as  an 
analysis  tool  will  be  investigated  through  use  of  both  real-world  and  simulated  data. 
Then,  in  comparing  the  plotted  IMFs  resulting  from  employment  of  the  CEMD  method 
both  before  and  after  the  application  of  the  Hilbert  transform  using  various  mathematical 
transforms,  the  final  decision  will  be  made  about  the  use  of  the  CEMD  method  and  the 
HHT  as  analysis  tools. 

1.4.  Approach/Methodology 

The  methodology  that  will  be  used  to  complete  the  research  will  be  described  in 
greater  detail  in  the  methodology  chapter,  but  a  brief  overview  is  provided  in  the 
following  paragraphs. 

1.4.1.  Data  Collection 

The  first  data  set  evaluated  was  the  real-world  data  set  provided  by  the  sponsor. 
In  addition  to  this  data  set,  four  sets  of  simulated  data  collected  at  AFIT’s  RCS  range 
were  used  in  comparing  with  the  decomposed  real-world  data  set  after  the  CEMD  method 
is  employed.  These  four  sets  of  simulated  data  included:  a  cylinder  with  two  end  caps;  a 
cylinder  with  one  end  cap  and  one  open  end  (or  cavity  cylinder);  an  ogive;  and  a  dihedral 
comer  reflector.  Once  all  data  was  collected,  it  will  be  calibrated  using  the  MATLAB® 
graphical  user  interface  (GUI)  called  ALPINE©  [18].  Then,  it  will  be  converted  into  a 
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MAT-file  that  can  be  read  into  MATLAB®.  At  this  point,  the  algorithm  introduced  in 
the  following  paragraph  will  be  employed. 

1.4.2.  Explanation  of  Chosen  Algorithm 

Research  has  been  done  concerning  methods  that  presently  exist  in  the  realm  of 
analyzing  complex- valued  data  sets  that  extend  the  EMD  technique.  Of  the  three 
methods  that  are  most  recognized  by  the  academic  community,  the  method  that  was 
employed  in  the  data  analysis  for  this  work  was  the  BEMD  method  developed  by  Rilling 
and  Flandrin  [11].  The  BEMD  method  was  used  as  a  model,  but  there  were  some 
modifications  to  use  the  method  on  the  dynamic  object  data  set  provided.  Once  the 
BEMD  algorithm  was  created,  the  analysis  of  the  data  was  accomplished. 

After  the  BEMD  algorithm  was  refined  and  tested  on  code  provided,  the  real- 
world  and  simulated  data  were  analyzed  by  the  BEMD  code  developed  for  this  research 
effort.  After  the  data  was  analyzed  by  the  BEMD  technique  code,  three  different 
transforms  were  applied  to  the  resulting  IMFs.  First,  the  Hilbert  transform  was  applied  to 
the  decomposed  data.  Second,  the  FFT  was  applied  to  the  decomposed  data  both  before 
and  after  the  Hilbert  transform  application.  Finally,  the  windowed  FFT  was  applied  to 
the  decomposed  data  both  before  and  after  the  Hilbert  transform,  represented  as  a 
Doppler-Time-Intensity  (DTI)  plot.  For  all  three  applied  transforms,  the  IMFs  of  the 
complex  RCS  data  were  plotted  and  recorded,  being  used  in  the  data  analysis  portion  of 
the  research. 

1.4.3.  Qualitative  Data  to  be  Analyzed 

The  analyzed  data  was  qualitative  data,  in  the  form  of  graphs  created  in 
MATFAB®.  The  plotted  IMFs  resulting  from  the  application  of  the  Hilbert  transform, 
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the  FFT,  and  the  Windowed  FFT  were  the  graphs  compared  for  data  analysis.  These 
plotted  IMFs  were  compared  and  distinct  features,  to  include  differences,  similarities,  and 
unique  points,  were  documented.  The  sets  of  plotted  IMFs  from  the  three  mathematical 
transformed  mentioned  earlier  in  this  paragraph  were  the  data  analyzed. 

1.5.  Overview  of  Chapters 

The  remaining  chapters  of  the  thesis  include  the  following  components:  the 
literature  review  section;  the  methodology  section;  the  data  analysis  and  results  section; 
and  the  conclusions  section.  The  literature  review  will  contain  analysis  of  work  already 
done  on  the  topic  of  the  EMD  technique  and  the  HHT,  including  articles  relevant  to  the 
research  for  the  problem  being  investigated.  The  methodology  section  details  the 
approach  to  complete  the  research  and  method  employed  to  create  the  qualitative  data  to 
be  analyzed  in  the  data  analysis  section.  The  main  idea  of  the  methodology  is  to 
qualitatively  represent  the  problem  as  graphs,  as  well  as  provide  a  way  to  assess  the  use 
of  the  CEMD  method  and  the  HHT  as  analysis  tools.  After  the  methodology  is 
explained,  the  plotted  data  sets  will  be  discussed  and  analyzed  in  the  data  analysis 
section.  The  resulting  plotted  IMFs  created  by  the  data  will  be  analyzed  and  compared  in 
order  to  determine  whether  the  CEMD  algorithm  and  the  HHT  provide  enhanced  fidelity 
than  the  FFT.  The  thesis  will  continue  with  a  discussion  of  the  data  created  by  the 
CEMD  coding.  Concluding  thoughts  will  be  made  in  the  final  section  as  to  the 
assessment  of  the  CEMD  algorithm  and  the  HHT  as  possible  analysis  tools,  as  well  as 
how  the  Hilbert  transform  compares  with  the  FFT  analysis  of  dynamic  object  data  sets. 


11 


II.  Literature  Review 


2.1.  Introduction 

Traditionally,  the  use  of  “Fourier  spectral  analysis”  [1]  through  application  of  the 
Fourier  transform,  as  well  as  the  faster  algorithm  of  the  FFT,  has  been  used  to  quantify 
various  sets  of  data  in  the  signal  analysis  department.  The  restrictive  nature  of  the 
Fourier  transform  primary  need  for  a  linear  system  and  strictly  periodic  or  stationary  data 
[1],  other  transforms  have  been  researched  as  possible  supplemental  or  replacement 
analysis  tools  in  the  fields  where  nonlinear  and  non-stationary  data  exist. 

More  recently,  investigation  into  a  newer  mathematical  transforms  have  been 
conducted.  One  transform  called  the  HHT  was  introduced  by  Dr.  Huang  et  al.  [1]  and 
has  since  been  used  in  numerous  fields  of  science  and  engineering,  ranging  from  seismic 
analysis  to  heartbeat  patterns  [3],  [4],  [6],  and  [14].  The  HHT  consists  of  two 
components:  (1)  the  Empirical  Mode  Decomposition  (EMD)  method  and  (2)  the  Hilbert 
transform. 

The  EMD  method  is  an  adaptive  technique  in  which  a  given  signal  is  decomposed 
into  a  set  of  oscillating  components  through  use  of  the  sifting  process.  The  oscillatory 
components  are  called  the  IMFs  of  the  given  signal  and  are  representations  of  the 
oscillating  nature  embedded  in  the  data.  The  FFT  has  a  defined  set  of  basis  functions; 
however,  the  IMFs  that  are  decomposed  from  the  data  are  considered  the  basis  functions 
for  the  signal.  The  EMD  technique  and  application  of  the  Hilbert  transform  to  create  the 
Hilbert  spectrum  comprise  the  HHT. 

The  EMD  accepts  only  real-valued  signals  and  could  lead  to  the  loss  of  some 
information  provided  by  the  signal.  Therefore,  research  has  been  performed  in  the  past 
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five  years  to  create  a  way  to  extend  the  EMD  technique  to  analyze  complex-valued 
signals,  consisting  of  the  real-valued  and  imaginary-valued  data.  Several  researchers 
[9,10,11]  have  also  developed  methods  that  fall  under  the  topic  of  the  CEMD;  they  have 
extended  the  EMD  technique  to  incorporate  complex-valued  data  sets  rather  than  only  the 
real-valued  components  of  the  data. 

2.2.  Literature  Review  Structure 

The  work  accomplished  by  Dr.  Huang  et  al.  [1]  is  analyzed  by  the  various  ways 
the  HHT  and  EMD  technique  have  been  adapted  into  a  number  of  fields  of  study  [3] -[4], 
[6],  [14]  and  will  conclude  with  the  more  recent  developments  for  the  CEMD  technique 
[9]-[l  1],  [15]-[17],  as  well  as  multivariate  signals  [12].  Over  the  course  of  approximately 
15  years  of  research,  many  breakthroughs  and  discoveries  concerning  the  use  of  the  HHT 
have  occurred,  more  specifically  concerning  the  EMD  technique  inherent  to  the  HHT. 
Upon  completing  this  literature  review,  a  greater  understanding  of  the  HHT  and  EMD 
technique  is  expected. 

2.3.  The  Hilbert-Huang  Transform 

Where  the  details  of  the  HHT  are  concerned,  Dr.  N.  E.  Huang  is  a  leading 
authority  in  the  field.  Huang  et  al.  wrote  introduced  “[t]he  empirical  mode 
decomposition  and  the  Hilbert  spectrum  for  nonlinear  and  non- stationary  time  series 
analysis”  [1]  in  1998  and  introduced  a  new  method  for  analyzing  nonlinear  and  non¬ 
stationary  data,  known  as  the  HHT.  There  are  two  components  of  the  HHT,  as  stated 
above,  but  the  key  part  of  the  HHT  is  the  EMD  technique.  Huang  et  al.  claimed  that, 
through  the  application  of  the  EMD  technique,  “any  complicated  data  set  can  be 
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decomposed  into  a  finite  and... small  number  of  ‘intrinsic  mode  functions’  that  admit 
well-behaved  Hilbert  transforms”  [1].  Much  interest  has  been  shown  for  use  of  the  EMD 
technique  due  to  its  adaptive  nature  and  the  claim  that  it  is  also  “highly-efficient”  [1]  due 
to  that  adaptiveness.  Another  important  trait  the  EMD  technique  possesses  is  its  ability 
to  be  applied  to  a  localized  region.  In  order  to  apply  the  EMD  to  nonlinear  and  non¬ 
stationary  time  series,  the  necessary  conditions  of  the  data  are  that  of  (1)  locality  and  (2) 
adaptivity. 

2.3.1.  Pre-Existing  Non- Stationary  Methods 

In  preparing  for  the  explanation  of  the  development  of  the  EMD  technique  and 
Hilbert  spectrum,  Huang  et  al.  present  various  pre-existing  non- stationary  data  processing 
methods  [1].  The  following  methods  work  for  non-stationary  (or  non-periodic)  data,  but 
depend  heavily  on  Fourier  analysis,  leading  the  applications  of  these  methods  to  be 
limited  to  linear  data. 

The  first  method  discussed  is  the  most  basic  method  called  the  spectrogram. 
Huang  et  al.  claim  that  the  spectrogram  “is  nothing  but  a  limited  time  window-width 
Fourier  spectral  analysis”  [1].  Because  it  depends  on  Fourier  analysis,  it  is  not  used  for 
the  analysis  of  nonlinear  and  non-stationary  data. 

The  next  approach  is  the  wavelet  analysis,  an  adjustable  window  Fourier  analysis, 
defined  by  the  following  equation: 


W(a,b;X,y/)  =  \a\  [  X{t)y/ 

'  1  J  -oo 


^  t-b^ 


(5) 


dt 


V  a  J 

where  !//*(•)  is  the  wavelet  function,  a  is  the  dilation  factor  and  b  is  the  translation  of 


the  origin.  The  physical  explanation  of  (5)  is  that  W(a,b\X, y/')  is  the  “energy”  of  X  of 
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scale  a  at  t  =  b[  1].  The  continuous  wavelet  analysis  is  of  an  analytic  form,  in  which  the 
problem  employing  it  occurs  with  the  Morlet  wavelet  as  an  example,  making  the 
“quantitative  definition  of  the  energy- frequency-time  distribution  difficult”  [1].  Even 
though  the  wavelet  analysis  contains  difficulties,  Huang  et  al.  use  the  wavelet  analysis  in 
the  validation  of  the  Hilbert  spectrum. 

A  third  method,  called  the  Wigner-Ville  distribution,  is  also  sometimes  referred  to 
as  the  Heisenberg  wavelet.  By  definition,  the  Wigner-Ville  distribution  is  the  Fourier 
transform  of  the  central  covariance  function,  defined  by  the  following  equation: 

<•00 

V(cD,t)=\  Cc(j,t)e~mTdr  (6) 

J  —00 

where 


(  l  }  1  "j 

C(r,t)  =  X  t  —  r  X  t  +  —  r 
V  2  )  ^  2 


(7) 


The  problem  with  the  Wigner-Ville  distribution  occurs  with  the  cross  terms  that 
result  from  the  negative  energy  components,  thus  resulting  in  a  windowed  Fourier 
analysis  [1].  The  Fourier  analysis  limitations  are  forced  upon  the  Wigner-Ville 
distribution  analysis,  primarily  the  linearity  condition. 

Another  method  is  the  evolutionary  spectrum,  where  the  classical  Fourier  analysis 
is  extended  to  a  more  generalized  basis.  Thus,  a  method  is  sought  to  define  the  basis, 
{(j){co,t)] ,  without  defining  it  prior  to  the  application  of  the  method. 


The  final  method  introduced  by  Huang  et  al.  prior  to  the  EMD  technique  and 
Hilbert  spectrum  was  the  empirical  orthogonal  function  (EOF)  expansion.  EOF 
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expansion  states,  for  any  real  signal,  z(x,t) ,  the  application  of  the  EOF  will  reduce  the 


signal  to: 


z(x,t)  =  YJl=iak(t)fk(x)  (8) 

The  expansion  basis  for  EOF  is  derived  from  the  data,  showing  it  as  adaptive. 
Notwithstanding,  its  main  problem  is  the  uncertainty  of  its  true  meaning  where  non- 
stationarity  and  nonlinearity  are  concerned,  alluding  to  the  conclusion.  EOF  is  not  an 
effective  improvement  from  those  methods  dependent  upon  Fourier  analysis. 

2.3.2.  Method  Developed  by  Huang  et  al. 

The  method  introduced  by  Huang  et  al.  is  a  general  method  that  consists  of  two 
steps.  The  first  is  the  decomposition  of  the  data  through  the  EMD  technique,  where  the 
data  are  decomposed  into  a  number  of  IMFs.  The  next  step  is  to  apply  the  Hilbert 
transform  to  those  IMFs,  resulting  in  the  Hilbert  spectrum.  Contained  within  the  Hilbert 
spectrum  are  the  instantaneous  (or  local)  frequency  and  instantaneous  energy,  and  are 
used  for  analysis  rather  than  the  global  frequency  and  energy  as  defined  by  Fourier 
analysis  [1]. 

2.3.3.  Intrinsic  Mode  Functions 

Huang  et  al.  proposed  a  “class  of  functions  designated  as  intrinsic  mode 

functions,”  [1]  where  the  formal  definition  is  as  follows: 

“An  intrinsic  mode  function  (IMF)  is  a  function  that  satisfies  two  conditions:  (1) 
in  the  whole  data  set,  the  number  of  extrema  and  the  number  of  zero  crossing 
must  either  equal  or  differ  at  most  by  one;  and  (2)  at  any  point,  the  mean  value  of 
the  envelope  defined  by  the  local  maxima  and  the  envelope  defined  by  the  local 
minima  is  zero”  [1]. 
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The  first  condition,  as  stated  above,  is  similar  to  the  narrowband  requirements  for 
a  stationary  process.  The  second  condition  presented  takes  the  global  requirement  from  a 
stationary  data  set  and  adapts  it  to  the  local  level  for  non-stationary  data  sets.  The  ideal 
requirement  is  for  the  local  mean  of  the  given  data  set  to  be  zero.  The  IMFs  represent  the 
oscillatory  modes  that  are  embedded  in  the  data  and  can  be  in  the  form  of  either 
amplitude-modulated  signals  or  frequency-modulated  signals. 

2.3.4.  The  Empirical  Mode  Decomposition  Method 

The  EMD  method  is  also  commonly  referred  to  as  the  sifting  process.  This  new 
method  is  used  with  both  nonlinear  and  non-stationary  data.  In  order  to  help  with  the 
analysis  of  nonlinear  and  non-stationary  signals,  the  EMD  method  has  successfully  been 
used  in  applications  of  various  disciplines  due  to  its  versatile  data-driven  signal  analysis 
ability.  The  EMD  method  is  a  new  technique,  pioneered  specifically  for  the  purposes  of 
adaptively  representing  nonlinear  and  non-stationary  signals  as  sums  of  zero-mean 
amplitude-  and  frequency-modulated  components.  The  one  new  feature  of  this  method, 
when  compared  to  previously  existing  methods,  is  that  the  EMD  technique  is  “intuitive, 
direct,  a  posteriori  and  adaptive,  with  the  basis  of  the  decomposition  based  on... the 
[given]  data”  [1], 

The  main  idea  of  the  EMD  technique  is  to  decompose  the  signal  into  its 
oscillatory  modes.  As  long  as  the  two  conditions  stated  in  the  IMF  definition  in  section 
2.2.3  are  satisfied,  the  EMD  method  can  use  the  envelopes  defined  by  the  maxima  and 
minima  separately.  After  the  extrema  are  identified,  the  local  maxima  are  connected  with 
a  cubic  spline  as  the  upper  envelope.  The  same  process  is  performed  for  the  lower 
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envelope.  Sifting  is  a  similar  process,  in  which  the  finest  local  mode  is  separated  from 
the  rest  of  the  data  [1], 

The  EMD  method  is  designed  to  “reduce  non-stationary,  multicomponent  signals 
to  a  series  of  amplitude-  and  frequency-modulated  contributions”  [6]  and  can  be  used  to 
gain  significant  information  inherent  to  the  signal.  Although  other  methods  exist  for  non- 
stationary  analysis,  the  EMD  method  differs  from  wavelet  decomposition  in  which  the 
“filters  of  the  filter  band  do  not  correspond  to  sub-band  filtering  but  instead  to  signal- 
dependent,  time-variant  filters”  [6]. 

The  EMD  technique  operates  in  the  time-domain  and  adaptively  decomposes  a 
signal  into  a  set  of  basis  functions  called  the  IMFs,  and  data  can  be  considered  to  be 
mapped  onto  a  space  spanned  by  the  IMFs  [10].  By  applying  the  Hilbert  transform  to  the 
IMFs,  the  “instantaneous  frequency”  is  introduced. 

To  apply  the  technique  of  EMD,  the  given  signal  must  be  considered  at  the  local 
oscillation  level.  The  algorithm  that  Flandrin  et  al.  employs  is  the  algorithm  described  in 
Table  2  [8], 


Table  2:  Flandrin  et  al.  EMD  Algorithm 


Given  a  signal  x(t) : 

1. 

Identify  all  extrema  of  x(t) 

2. 

Interpolate  between  minima  and  respective  maxima,  ending  up  with  some 
envelope  emm  (t)  and  its  respective  emax  (1 ) . 

3. 

Compute  the  mean:  m(t)  =  ^  (emm  ( t )  +  emax  (t)) 

4. 

Extract  the  detail:  d(t)  =  x(t)-m(r) 

5. 

Iterate  on  the  residual:  m(t) 
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In  the  above  algorithm,  the  detail,  denoted  by  d(t),  corresponds  to  the  oscillation 


terminating  at  the  two  minima  and  passing  through  the  maximum  value  of  the  oscillating 
wave,  existing  between  the  two  extrema.  The  local  trend,  denoted  by  m(t)  ,  corresponds 
to  the  low-frequency  part.  The  sifting  process  is  employed  on  the  signal,  where  an  IMF 
and  residual  are  extracted  and  the  iterative  algorithm  is  performed  on  the  residual  portion. 
The  sifting  process  is  applied  to  d(t )  until  dir)  can  be  considered  as  zero-mean, 
according  to  the  stopping  criterion  [8], 

The  method  of  EMD  considers  a  signal  at  the  scale  of  its  local  oscillations  and 
attempts  to  formalize  the  idea  that  “signal  =  fast  oscillations  superimposed  on  slow 
oscillations”  [11].  The  EMD  technique  is  designed  to  define  local  “low  frequency” 
components  as  the  local  trend,  n\[x\(t) ,  where  this  local  trend  then  supports  a  local  “high 

frequency”  component  as  a  zero-mean  oscillation,  or  local  detail,  d^x\(t)  . 

The  signal  is  represented  by  the  following  expression: 

x(t)  =  mx[x](t)  +  dx[x](t)  (9) 

where  dJxKO  corresponds  to  an  IMF.  The  sifting  process  is  performed  on  (9)  and  the 
signal  expression  becomes: 

40  =  mk  [x](0  +  dk  [xj(0  (10) 

Once  the  convergence  criterion  has  been  met,  the  local  detail  and  local  trend  are 
represented  as,  dx[x\(t)  -  S"[x](0  and  mx[x\(t)  =  x(t)~ d\x](t),  respectively  [11]. 
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2.3.5.  The  Sifting  Process 

Huang  et  al.  suggested  a  data-adapted  method  in  which  an  oscillating  wave  is 
extracted  from  a  given  signal  [13].  Each  oscillating  wave  is  defined  as  an  IMF  that 
satisfies  the  two  conditions  outlined  by  Huang  et  al.  [1].  The  sifting  process  is  an 
iterative  process,  meaning  that  if  the  two  conditions  outlined  by  Huang  et  al.  are  not  met, 
then  the  sifting  procedure  will  be  repeated  until  the  two  conditions  are  satisfied  [13].  A 
“stopping  rule”  [13]  is  applied  when  all  that  remains  of  the  original  input  signal  is  the 
residual  after  the  IMFs  have  been  extracted  from  the  given  signal. 

2.3.6.  The  Hilbert  Transform 

For  a  real  signal,  x(t) ,  the  analytic  signal  is  defined  as 

z(t)  =  x(t)  +  jy(t)  (11) 

In  (11),  y(t)  represents  the  Hilbert  transform  of  the  real  signal,  where 


y(t)  =  — W  ^-ds  (12) 

n  J~™t-s 

and  P  is  the  Cauchy  principal  value.  To  describe  the  instantaneous  frequency  in  its 
correct  form,  the  analytic  signal  must  be  defined  by  polar  coordinates,  leading  to  the 
analytic  signal  being  defined  as: 

z(t)  =  a(t)exp(j0(t))  (13) 

In  (13),  a(t)  is  the  amplitude  and  defined  as  follows: 
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and  6(t)  represents  the  phase,  defined  by 


0(t)  -  tan 


ryi 

x(t)y 


(15) 


Finally,  the  instantaneous  frequency  as  time-varying  phase  is  defined  as, 


dO(t) 

clt 


(16) 


2.3. 7.  The  Hilbert  Spectrum 

When  given  a  non-stationary  signal  with  variable  frequency  and  amplitude  change 
over  a  period  of  time,  there  is  a  need  to  have  a  more  adaptive  and  flexible  notion  of 
frequency.  The  concept  of  the  Hilbert  spectrum  and  instantaneous  frequency  were 
detailed  by  Huang  et  al.  [1]  through  the  application  of  the  Hilbert  transform. 

Once  the  IMFs  have  been  decomposed  from  the  original  signal,  the  Hilbert 
transform  can  be  applied  to  each  of  the  IMFs  individually,  resulting  in  the  components  of 
the  Hilbert  spectrum.  After  applying  the  Hilbert  transform  and  computing  the 
instantaneous  frequency  using  (16),  the  data  set  can  now  be  expressed  as: 


x  (o = XL  ak (*)  exp  ( /j  ojk  (odt)  (i7) 

where  ak  represents  the  amplitude  of  each  component  as  a  function  of  time  and  cok 
represents  the  instantaneous  frequency  of  each  component  as  a  function  of  time. 
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In  comparison,  the  Fourier  representation  would  be  expressed  as: 


=  (18) 
where  the  amplitude  ak  and  frequency  cok  are  both  constants,  rather  than  variable  with 
time. 

By  using  the  IMF  expansion  rather  than  the  Fourier  expansion,  the  restrictions  of 
expansion  on  a  linear  and  stationary  data  set  disappears  and  the  function  can  now  handle 
variable  amplitude-  and  frequency-modulation,  leading  to  an  easier  analysis  of  the 
nonlinear  and  non- stationary  data  sets. 

Once  the  instantaneous  frequency  is  calculated,  the  Hilbert  spectrum  can  be 
represented  by  the  triplet  of  {r,&>.(r), 4(0}  (°r  time-frequency-amplitude)  in  the  time- 
frequency  plane,  where  A;  (7)  is  the  amplitude  of  the  analytic  signal,  corresponding  to  the 
instantaneous  frequency  cof(t) . 

The  original  EMD  performs  the  mapping  expressed  below: 

x\n\  =  k[n\  +  r[n\  (19) 

From  this  expression,  the  IMFs,  denoted  by  dk[n\,  represent  a  unique  time- 
frequency  analyzer  allowing  for  analysis  of  the  instantaneous  frequency.  The 
combination  of  the  concept  of  instantaneous  frequency  and  the  EMD  technique  makes  the 
EMD  framework  so  powerful  for  time-frequency  signal  analysis  [9]. 
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2.4.  Algorithm  Variations 

In  the  article  “On  Empirical  Mode  Decomposition  and  its  Algorithms,”  [5]  by 
Rilling  et  al.,  a  new  data-driven  technique  of  EMD  is  presented  and  issues  related  to  its 
effective  implementation  are  discussed.  The  technique  of  EMD  is  faced  with  the 
difficulty  of  not  having  an  analytic  form,  being  defined  only  by  an  algorithm,  making 
theoretical  analysis  and  performance  evaluation  nearly  impossible  [5],  In  addition  to 
presenting  the  problems  inherent  in  the  EMD  technique,  Rilling  et  al.  also  propose  some 
variations  on  the  EMD  algorithm.  Results  were  obtained  from  numerical  simulations  in 
order  to  support  an  interpretation  of  the  method  in  terms  of  adaptive  constant-Q  filter 
banks  [5]. 

2.4.1.  Algorithmic  Variations 

The  aim  of  Rilling  et  al.  in  presenting  these  algorithmic  variations  [5]  was  to 
make  the  choices  made  by  the  user  more  precise  and  to  recommend  specific  rationales 
behind  the  decisions  that  the  user  makes  prior  to  implementing  the  EMD  algorithm. 
There  are  two  variations  that  Rilling  et  al.  present  in  this  paper.  The  first  is  the  Local 
EMD  and  the  second  is  the  On-line  EMD. 

The  Local  EMD  algorithmic  variation  includes  a  variation  made  concerning  the 
initial  EMD  algorithm  formulation.  The  authors  introduce  an  intermediate  step  in  the 
sifting  process,  where  the  large  error  zeros  exist.  The  sifting  process  eliminates  the 
problem  of  over-iterating  the  entire  signal  by  targeting  the  zeros  that  cause  the  largest 
error  [5].  In  the  algorithm,  the  extra  iterations  performed  are  denoted  by  a  weighting 
function,  w(t). 
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This  weighting  function  is  introduced  in  Step  4  of  the  original  EMD  algorithm 
detailed  in  Table  2,  such  that  the  new  definition  of  the  detail  is  now: 

d (t)  =  x(t)  -  w(t )m(t )  (20) 

The  second  variation,  referred  to  as  the  On-line  EMD,  is  based  on  the  fact  that  the 
sifting  step  relies  on  the  interpolation  between  the  local  extrema.  Rilling  et  al.  claim  that 
because  the  interpolation  of  the  extrema  requires  only  a  finite  number  of  interpolations, 
“that  the  extraction  of  a  mode  could  therefore  be  possible  blockwise,  without  the 
necessary  knowledge  of  the  whole  signal  (or  previous  residual)”  [5].  In  order  to  realize 
this  algorithmic  variation,  a  sliding  window  was  implemented  on  the  original  EMD 
algorithm,  as  described  in  Table  4. 

2.4.2.  Performance  Elements 

In  addition  to  introducing  the  previous  two  algorithmic  variations  on  the  EMD 
algorithm.  Rilling  et  al.  also  identify  some  performance  elements  causing  problems  when 
implementing  the  EMD  algorithm.  Because  the  EMD  technique  is  defined  by  an 
algorithm,  performance  evaluation  is  difficult  and  requires  simulation  experiments  [5]. 
The  first  problem  they  address  is  the  idea  of  tones  and  sampling.  The  EMD  is  expected 
to  be  the  identity  operator  with  only  one  tone  and  no  residual;  however,  in  actuality,  this 
is  not  a  true  statement.  The  issue  arises  from  the  fact  that  tone  estimation  depends 
heavily  on  the  tone  frequency,  the  application  of  the  EMD  technique  results  in  a  number 
of  IMFs,  as  well  as  a  residual  component. 
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2.5.  Empirical  Mode  Decomposition  as  a  Filter  Bank 

Flandrin  et  al.  presented  [8],  in  which  they  reported  on  experiments  involving 
fractional  Gaussian  noise  to  better  understand  how  the  EMD  technique  behaves  in 
situations  involving  broadband  noise.  They  conclude  that  the  technique  of  EMD  acts  as  a 
dyadic  filter  bank,  resembling  filters  existing  in  wavelet  decomposition  [8]. 

2.5.1.  Fractional  Gaussian  Noise 

The  final  form  of  the  EMD  results  in  the  representation  as  follows: 

x(t)  =  YJ^]c[(t)  +  r(t)  (21) 

where  r(t )  stands  for  the  residual  trend  of  the  entire  signal  and  ct(t)  represents  the  IMFs 
throughout  the  decomposition  [8]. 

Less  attention  has  been  given  to  realistic  situations  involving  noise,  where  most 
studies  have  been  performed  based  on  simulations,  resulting  in  less  understanding  of  the 
decomposition  that  EMD  can  achieve  when  applied  to  a  stochastic  (or  intrinsically  non- 
deterministic)  process.  In  order  to  help  explain  applying  the  EMD  technique  to  such 
processes  that  involve  noise,  Flandrin  et  al.  performed  extensive  realistic  simulations  in 
order  to  show  that  the  EMD  technique  performs  like  a  dyadic  filter  when  applied  to  noise 
processes  [8]. 

2.5.2.  Flandrin  et  al.  Conclusions 

As  a  result  of  their  research  efforts,  Flandrin  et  al.  were  able  to  report  on  the  “first 
numerical  experiments  aimed  at  supporting  the  claim:  “. .  .the  built-in  adaptivity  of  EMD 
makes  it  behaves  spontaneously  as  a  ‘wavelet-like’  filter  bank”  [8].  The  technique  of 
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EMD  naturally  copes  when  there  are  superimposed  IMFs  and  the  benefits  of  the  EMD 
method  are  similar  to  those  of  wavelet-based  methods  [8]. 

2.6.  EMD  Applications 

The  applications  described  in  the  following  section  employ  only  real-valued  data, 
rather  than  complex-valued  data. 

2.6.1.  Seismic  Application 

Magrin-Chagnolleau  and  Baraniuk  [7]  propose  a  new  technique  called  the  EMD 
is  described  and  applied  to  the  investigation  of  a  seismic  trace,  where  the  IMFs  and 
instantaneous  frequency  were  studied.  They  also  applied  the  EMD  technique  to  a  seismic 
section,  resulting  in  new  time-frequency  attributes. 

The  topic  in  which  Magrin-Chagnolleau  and  Baraniuk  applied  the  EMD  technique 
was  that  of  seismic  signals.  Much  like  many  real-world  signals,  seismic  signals  have  the 
property  that  they  are  non-stationary  [7].  Due  to  the  inability  of  Fourier  analysis  to 
analyze  non-stationary  and  nonlinear  signals,  Fourier  analysis  provides  unsatisfying 
results  due  to  the  frequency  changes  that  occur  with  respect  over  time  of  the  seismic 
signal.  Because  the  EMD  method  is  an  adaptive  decomposition  technique  that 
decomposes  the  signal  into  its  oscillating  components,  this  new  technique  has  potential  in 
analyzing  seismic  signals  when  compared  with  Fourier-based  analysis  tools. 

In  their  paper,  Magrin-Chagnolleau  and  Baraniuk  proposed  a  new  way  of 
decomposing  a  seismic  trace  into  its  IMFs  and  extracting  the  instantaneous  frequency  of 
each  IMF.  The  next  step  in  their  research  would  be  to  extract  other  time-frequency 
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attributes  based  on  different  calculations  in  the  time-frequency  plane  represented  by  the 
triplet,  4(0}  . 

2.6.2.  Heart  Rate  Variability  Application 

Balocchi  et  al.,  in  “Deriving  the  respiratory  sinus  arrhythmia  from  the  heartbeat 
time  series  using  Empirical  Mode  Decomposition,”  explored  an  application  of  the  EMD 
technique  [3].  Heart-rate  variability  (HRV)  is  a  well-known  phenomenon  and  is  of  great 
clinical  relevance  in  pathophysiologic  investigations.  However,  analyzing  HRV  is 
difficult  because  it  is  the  result  of  many  nonlinear  interacting  processes.  Any  linear 
analysis  tool  that  is  applied  to  the  HRV  has  the  potential  of  underestimating  or  missing 
information.  Therefore,  researchers  have  applied  EMD  analysis  to  decompose  the 
heartbeat  interval  series  into  their  IMFs  in  order  to  identify  the  modes  associated  with 
breathing  [3].  For  comparison  purposes,  Balocchi  et  al.  recorded  the  respiratory  signal 
simultaneously  with  the  tachogram  (or  EKG)  signal. 

As  previously  stated,  the  EMD  method  allows  the  analysis  of  nonlinear  and  non¬ 
stationary  time  series  through  the  analysis  of  their  IMFs.  In  this  application  performed 
by  Balocchi  et  al.,  the  authors  were  able  to  demonstrate  the  association  of  the  first  IMF 
extracted  from  a  tachogram  with  the  simultaneously  recorded  respiratory  signal  [3]. 

2.6.3.  Seismic  Reflection  Application 

The  article  “Application  of  the  Empirical  Mode  Decomposition  and  Hilbert- 
Huang  Transform  to  seismic  reflection  data,”  written  by  Battista  et  al.,  applied  the 
technique  of  EMD  to  the  study  of  seismic  reflection  data  [6].  There  have  been 
advancements  in  the  field  of  signal  processing  providing  possibly  improved  imaging  and 
analysis  of  “complex  geologic  targets  found  in  seismic  reflection  data”  [6],  The  EMD 
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technique  has  yet  to  be  recognized  as  a  standard  analysis  tool  by  the  seismic  community. 
Therefore,  the  reasoning  behind  this  experiment  is  to  demonstrate  the  ability  of  the  EMD 
technique  and  HHT  to  improve  seismic  reflection  data  quality. 

The  HHT  allows  for  signals  that  are  described  as  stochastic  (or  intrinsically  non- 
deterministic)  processes  to  be  analyzed  by  using  instantaneous  attributes,  such  as 
frequency  or  displacement,  in  the  time-frequency  domain.  Two  reasons  the  authors 
applied  the  HHT  to  the  data  were:  (1)  to  assess  the  ability  of  the  EMD  and  HHT  to 
quantify  geologic  information  in  the  time  and  time-frequency  domain  and  (2)  to  develop 
superior  filters  by  using  the  instantaneous  attributes.  The  main  objective  of  the 
experiment  was  to  determine  whether  HHT  allows  for  filter  design  using  its  empirically- 
derived  attributes  [6] . 

For  this  application,  the  HHT  was  first  used  to  compare  the  filtering  in  time- 
frequency  domain  against  that  of  the  frequency  domain  using  Fourier  transform.  Then, 
the  instantaneous  attributes  of  the  HHT  were  compared  to  those  produced  by  the  Hilbert 
transform,  where  the  EMD  technique  was  not  performed  with  the  Hilbert  transform.  By 
performing  these  two  comparisons  with  well-known  transforms,  the  authors  were  able  to 
demonstrate  the  strength  of  using  the  time-frequency  domain  filtering  and  the  necessity 
of  using  EMD  with  the  Hilbert  transform  [6] . 

EMD  and  HHT  were  not  presented  as  a  replacement  for  existing  methods,  but  the 
objectives  of  the  study  were  met.  The  HHT  is  an  impressive  analysis  tool  due  to  its 
ability  to  preserve  phase  and  amplitude  while  empirically  separating  the  signal  from 
noise.  Battista  et  al.  determine  that  future  goals  include  integrating  HHT  with 
“amplitude-versus  offset  processing  of  gas  hydrates”  [6], 
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2.6.4.  Cycle  and  Trend  Mode  Application 

Ehlers  and  Way  discussed  the  usefulness  of  an  objective  scientific  approach  for 
the  identification  of  cycle  or  trend  modes  in  the  market  [14].  While  a  number  of  tools  are 
already  available  to  provide  distinction  between  the  two  modes,  in  this  paper  a  unique 
new  approach  will  be  used  to  help  determine  the  market  mode  called  Empirical  Mode 
Decomposition. 

Ehlers  and  Way  determined  that  cycle  mode  components  of  market  activity  can  be 
identified  using  a  band  pass  filter.  An  uptrend,  which  represents  a  cycle  market  mode, 
can  be  identified  as  the  positive  average  of  the  filtered  data  over  cycle  periods  and  in  a 
similar  manner,  a  downtrend,  representing  a  trend  market  mode,  is  identified  by  the 
negative  average  of  the  filtered  data  over  cycle  periods  [14].  Finally,  the  delineation 
between  cycle  and  trend  modes  can  be  made  by  the  trend  line  deriving  using  the 
Empirical  Mode  Decomposition. 

2.6.5.  Petrophysical  Model  Application 

Huang  and  Milkereit  explore  another  use  for  the  EMD  method  because  of  the 
importance  of  spectral  analysis  for  seismic  data  processing  and  interpretation  [4].  Due  to 
the  fact  that  the  frequency  contents  of  seismic  data  vary  with  time,  the  medium  is  a  non¬ 
stationary  one.  The  advantage  of  the  HHT  and  EMD  is  they  do  not  require  presumed  set 
of  functions  as  previous  methods,  allowing  the  projection  of  non-stationary  and  nonlinear 
signals  onto  a  time-frequency  plane  using  the  Intrinsic  Mode  Functions.  Huang  and 
Milkereit  compare  the  findings  of  applying  EMD  method  with  those  of  applying  the 
wavelet  transform  and  the  S-transform.  In  previous  papers  in  the  seismic  community, 
spectral  decomposition  in  seismic  exploration  produced  a  continuous  time-frequency 
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expression  of  a  seismic  trace  [4].  The  EMD  method  generates  necessary  adaptive  bases 
from  data. 

Huang  and  Milkereit  present  two  applications  for  the  EMD  technique  in  their 
paper.  The  first  application  presented  is  the  synthetic  time  series,  where  they  applied  the 
EMD  method  to  three  time  series  similar  to  applications  used  in  three  pre-existing  papers 
[4].  After  the  first  application,  they  concluded,  the  EMD  method  provides  superior 
results  to  complex  wavelet  transform  and  S-transform  in  terms  of  temporal  and  spectral 
resolution. 

The  authors  also  applied  the  EMD  method  to  decompose  well-log  data.  The 
instantaneous  power  spectrum  density  function  provided  the  “depth  varying  stochastic 
properties  which  can  be  used  to  simulate  a  time  series  of  heterogeneous  medium  at  every 
depth”  [4],  Finally,  the  authors  determined  that  a  combination  of  two-dimensional  slices 
yields  a  heterogeneous  three-dimensional  earth  model  adaptive  to  non- stationary  along 
the  borehole.  In  conclusion,  Huang  and  Milkereit  decided  that  the  EMD  method  is  a 
helpful  tool  when  analyzing  the  seismic  data  used. 

2.7.  Extending  EMD  into  the  Complex  Domain 

Applications  in  the  previous  sub-section  employ  only  real-valued  data,  as  opposed 
to  complex-valued  data.  In  the  following  section,  complex-valued  data  was  used. 
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2.7.1.  Complex  Traces  and  Instantaneous  Frequency 

Given  a  real  signal  x(t),  the  corresponding  analytic  signal  (or  complex  trace)  is 
expressed  as: 

X(t)  =  x(t)  +  jH{x(t)},  (22) 

where  H{x(t) }  denotes  the  signal  corresponding  to  x(t)  and  was  obtained  using  the 
Hilbert  transform,  where  P  is  the  Cauchy  principle  value  of  the  integral  [7],  and  is 
expressed  as  follows: 

H{x(t)}  =  -pf  ^-dr  (23) 

Another  way  that  the  analytic  signal  can  be  obtained  is  detailed  below  in  Table  3  [7]. 

Table  3:  Analytic  Expression  of  Given  Signal 

Given  the  real  signal  x(t) : 

1.  Take  the  Fourier  transform  of  x(t) 

2.  Zero  the  amplitude  for  the  negative  frequencies  and  double  the  amplitude  for  the 

positive  frequencies 

3.  Take  the  inverse  Fourier  transform  of  the  resulting  signal 

In  polar  coordinates,  the  analytic  signal  can  be  expressed  as: 

X(t)  =  A(t)em\  (24) 

where,  A(t)  is  the  instantaneous  amplitude  and  0(t)  is  the  instantaneous  phase. 
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To  determine  the  instantaneous  radial  frequency,  the  following  equation  is 


implemented: 


co(t) 


dO(t) 

clt 


(25) 


2.8.  Methods  for  Empirical  Mode  Decomposition  in  Complex  Domain 

2.8.1.  Complex  Empirical  Mode  Decomposition 

Tanaka  and  Mandic  were  the  first  to  propose  a  method  for  extending  the  EMD 
method  into  the  complex  domain  [9].  They  proposed  to  achieve  the  complex  extension 
for  the  EMD  using  a  filter  bank  interpretation  of  the  EMD  mapping  and  by  use  of  the 
positive  and  negative  frequency  components  of  the  Fourier  spectrum.  This  method  yield 
complex-valued  IMFs,  facilitating  the  extension  of  the  standard  EMD  to  the  complex 
domain. 

The  authors  claim  that  the  EMD  method  is  a  “novel  signal  analysis  tool,  whereby 
the  underlying  notion  of  instantaneous  frequency  provides  an  insight  into  the  time- 
frequency  signal  features”  [9].  Through  the  research  accomplished  by  Huang  et  al.  [1], 
the  EMD  technique  is  an  established  tool  for  analyzing  non-stationary  and  nonlinear  data. 
Yet,  the  EMD  method  was  developed  only  for  real-valued  data,  leading  to  difficulties  in 
analysis  where  complex-valued  data  structures  exist. 

While  Tanaka  and  Mandic  claim  a  “simple  way”  to  extend  EMD  to  the  complex 
domain  would  be  to  apply  the  EMD  technique  separately  to  the  real  and  imaginary  parts 
of  a  complex-valued  signal;  however,  the  mutual  information  from  a  complex  quantity  is 
lost  when  the  signal  is  split  into  two  quantities  (real  and  imaginary).  Instead,  the  authors 
introduce  as  their  proposed  CEMD  the  concept  of  complex  IMFs  that  act  as  a  dyadic 
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filter  and  operate  directly  in  the  complex  domain,  where  the  signal  is  divided  into  the 
positive  and  negative  frequencies  [9].  The  only  requirement  for  this  method  is  that  the 
CEMD  preserves  the  filter  on  the  average  value. 

To  derive  the  CEMD,  the  complex-valued  data  set  must  first  be  decomposed  into 
its  positive  and  negative  frequency  components.  In  preparing  for  the  explanation  of  the 
method  Tanaka  and  Mandic  proposed,  let  {x[n]}  represent  a  complex-valued  time 

sequence  and  X(ejt0 )  represent  the  discrete-time  Fourier  transform  of  x[n\  e  C .  There 
are  two  possibilities  for  obtaining  the  desired  real  time  sequence  form  x[n\ ,  where  x[n\ 
is  generally  not  analytic,  making  one  of  the  above  mentioned  possibilities  unusable. 

The  other  possibility,  where  x[n\  is  not  assumed  to  be  analytic,  is  the  method 
used  for  extending  the  original  EMD  method  into  the  complex  domain.  To  extract  the 
positive  and  negative  frequency  component  from  x[n\ ,  an  ideal  bandpass  filter,  denoted 

by  BP(eJC0 ) ,  is  applied  to  the  original  signal. 

The  ideal  bandpass  filter  used  for  this  analysis  is  expressed  as 


BP(ejC0) 


J  1,0  <(0<K  1 

(o, ~7T  <  CO  <  Oj 


By  applying  the  bandpass  filter,  two  analytic  signals  are  generated: 


X+(e’(0)  =  BP(ejm)X(ejm) 

X_(ej(0)  =  BP{ej(°)X*  (eja) 

where  X*(eJW )  represents  the  complex  conjugate  of  the  signal. 


(26) 


(27,  28) 
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Once  the  bandpass  filter  has  been  applied,  the  inverse  Fourier  transform  is  employed: 


x+[n]  =  Re{F-1[x+(^'y)]} 

(  ,  r  .  Tl  (29,30) 

9 

where  S.H  denotes  the  real-part  of  the  signal  and  F  1  [•]  represents  the  inverse  Fourier 

transform  of  the  signal.  The  IMFs  can  be  obtained  using  the  following  summations  of 
(27)  and  (28): 


*+M  =  2],=i  x;M  +  r+M 

*An\  =  v  +  rAn\ 

9 

The  reconstruction  of  the  decomposed  complex  signal  is  as  follows: 


(31,32) 


x[n]  =  [x+[n]  +  jH[x+[n]])  +  {x_[n]  +  jH[x_[n]])  (33) 

9 

where  //[•]  represents  the  Hilbert  transform  of  the  signal.  To  obtain  the  i  th  complex 
IMF  of  the  complex  process  x[n] ,  the  following  equation  is  employed: 


>’,[«]  = 


x\n]  +  jH  [x([n]] ,  i  =  1, ...,  N+ 
(x;[«]  +  jff[x;[n]])  ,i  =  -N  -I 


(34) 


Therefore,  the  final  algorithm  representation  of  the  proposed  Complex  EMD  method  is: 


tin]  =  ,„ocW  +  An]  (35) 

Tanaka  and  Mandic  concluded  that  the  CEMD  method  can  be  achieved  based  on 
some  inherent  properties  of  the  complex  signals.  In  addition,  the  authors  have  been  able 
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to  apply  the  standard  EMD  to  corresponding  analytic  components  of  complex-valued 
data  used  in  their  paper  [9]. 

2.8.2.  Rotation  Invariant  Complex  Empirical  Mode  Decomposition 

Umair  Bin  Altaf  et  al.  propose  a  new  method  for  extending  the  EMD  technique 
into  the  complex  domain  is  proposed.  In  contrast  to  a  previous  method  proposed  by 
Tanaka  and  Mandic  [9],  this  method  is  achieved  in  a  generic  way  so  that  the 
mathematical  development  mirrors  that  of  the  original  EMD  method  [10].  Through  this 
method,  the  IMFs  are  complex  by  design  and  shown  to  provide  consistent  framework  for 
handling  real  and  imaginary  data. 

Traditional  time-frequency  analysis  methods  are  based  on  “a  priori”  mapping 
from  time  to  frequency  domains,  where  that  mapping  is  defined  by  “basis  functions” 
[10].  However,  this  “a  priori”  mapping  poses  problems  for  nonlinear  and  non-stationary 
signals  that  have  time  varying  statistical  characteristics,  and  a  single  basis  function  fails 
due  to  its  limited  accountability  for  the  variations.  In  addition,  this  “a  priori”  mapping 
also  compromises  the  physical  significance  of  the  signal  analysis  of  nonlinear  and  non- 
stationary  signals  [10]. 

The  authors  proposed  a  new  way  to  decompose  a  complex  signal  using  the 
method  of  EMD  and  was  achieved  by  making  use  of  the  complex  spline.  Using  the 
complex  spline  makes  it  possible  to  carry  out  the  arithmetic  and  algebraic  operations  of 
the  algorithm  in  the  complex  domain,  leading  to  a  single  set  of  IMFs  contained  in  the 
complex  domain  [10]. 

Before  describing  their  proposed  method,  the  authors  describe  the  method  that 
Tanaka  and  Mandic  proposed  CEMD  [9].  Concerning  their  method,  Umair  Bin  Altaf  et 
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al.  claim  this  method  is  a  “natural  and  generic  way  to  extend  EMD  to  the  complex 
domain  would  be  to  operate  in  the  complex  domain”  [10].  The  same  steps  for  the 
original  EMD  method  are  followed,  but  are  carried  out  in  the  complex  domain  with 
various  modifications. 

The  issue  that  arises  is  the  definition  of  an  extrema  in  the  complex  domain  and  a 
method  for  determining  it.  After  describing  a  number  of  definitions,  the  definition  of  an 
extrema  that  is  used  for  the  method  is:  “a  locus  where  the  angle  of  the  first  derivative 
(first-order  differential  vector)  changes  its  sign”  [10].  The  authors  have  assumed  that 
each  local  maximum  is  followed  by  a  local  minimum  and  vice  versa,  and  in  order  to 
prove  this  assumption  true,  the  average  of  the  envelopes  is  used.  Envelopes  for  the 
complex  signal  can  be  estimated  as  spline  interpolations  of  the  local  maxima,  and 
minima,  where  the  average  can  be  computed.  This  complex  spline  is  obtained  by 
computing  the  real  and  imaginary  parts  separately  [10]. 

The  final  steps  of  the  complex  algorithm  are  performed  like  the  original  EMD 
algorithm.  The  claimed  advantage  of  the  proposed  CEMD  method  is  that  it  does  not  split 
the  signal  into  two  parts  and  has  the  potential  to  be  extended  into  higher-dimensional 
cases  easily  [10],  unlike  the  method  proposed  by  Tanaka  and  Mandic  [9]. 

Umair  Bin  Altaf  et  al.  conclude  that  the  analysis  of  real-world  complex-valued 
data  shows  that  the  proposed  method  provides  new  insights  into  time-frequency  analysis 
of  nonlinear  and  non-stationary  signals,  which  was  not  possible  before  [10]. 

2.8.3.  Bivariate  Empirical  Mode  Decomposition 

Rilling  et  al.  present  a  new  method  for  extension  of  the  original  EMD  method  to 
the  complex  domain  in  their  paper  [11].  Initially,  the  method  of  EMD  was  limited  to  the 
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analysis  of  real-valued  time  series;  therefore,  an  extension  to  analysis  of  complex- value 
time  series  is  proposed  that  is  designed  to  extract  zero-mean  rotating  components  from 
the  signal,  where  the  original  EMD  extracts  zero-mean  oscillating  components  [11]. 

The  basic  idea  underlying  the  proposed  BEMD  is  that  a  “bivariate  signal  =  fast 
rotations  superimposed  on  slower  rotations”  [11].  The  slowly  rotating  component  has  to 
be  defined  as  the  mean  of  some  “envelope,”  where  the  envelope  is  represented  by  a  three- 
dimensional  tube  tightly  enclosing  the  signal.  Given  a  set  of  points  on  the  tube,  there  are 
at  least  two  ways  to  define  the  mean  of  the  envelope:  (1)  define  the  mean  as  the 
barycenter  of  four  points,  each  having  unit  mass  and  (2)  define  the  mean  as  the 
intersection  of  two  lines,  one  being  halfway  between  the  two  horizontal  tangents  and  the 
other  being  halfway  between  the  two  vertical  tangents  [11].  Due  to  the  second  definition 
being  less  prone  to  errors,  it  is  a  more  preferred  method.  The  goal  for  the  bivariate 
interpolation  is  the  same  as  with  the  original  EMD:  smooth  interpolation  with  as  few 
“bumps”  as  possible,  calling  for  the  use  of  a  cubic  spline. 

Given  an  angle  direction  that  performs  uniform  sampling  around  the  unit  circle, 
the  bivariate  extensions  are  defined  by  the  EMD  algorithm,  only  with  new  sifting 
elementary  operators  defined  by  Sm  and  SB2  .w'hich  correspond  to  the  algorithms  in 
Tables  4  and  5  [11]. 
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Table  4:  Algorithm  1  for  EMD  Bivariate  Extension 


2kn 

1.  Given  an  angle  direction:  pk  = -  for  1  <  k  <  N  complete  steps  2  through  4 

N 

2.  Project  the  complex-valued  signal  x(t) on  direction  pk  : 8  Pk  =  Re\_e]Pk x(t)\ 

3.  Extract  the  locations  {tf  }  (time  corresponding  to  the  angle  direction  and  IMF) 
of  the  maxima  of  8  „ 

Pk 

4.  Interpolate  the  set  {if ,  x(tf ) }  to  obtain  the  envelope  curve  in  direction 

Pk  ■  epk  (0 

5.  Compute  the  mean  of  all  envelope  curves:  mp)  =  —  ^  epk  (t) 

6.  Subtract  the  mean  to  obtain  SBlx(t )  =  x(t )  -  m(t) 

Table  5:  Algorithm  2  for  EMD  Bivariate  Extension 

2kn 

1.  Given  an  angle  direction:  pk  = -  for  1  <  k  <  N  complete  steps  2  through  4 

N 

2.  Project  the  complex-valued  signal  x(t)  on  direction  pk  :  8 pj  =  Re[c',/',/<  x</  )  ] 

3.  Extract  the  locations  [tf  ,df]  of  the  maxima  of  c/h 

4.  Interpolate  the  set  {tf  ,e^Pkd\ }  to  obtain  the  partial  (  or  partial  differential 
equation  of  the)  envelope  curve  in  direction  pk  :  e' p  (t) 

5.  Compute  the  mean  of  all  tangents:  m(t)  =  —  ^k  e'Pk  (t) 

6.  Subtract  the  mean  to  obtain  SB2x(t )  =  x(t)  -  m(t) 

The  BEMD  was  designed  so  that  the  signals  rotating  around  zero  are  the  outputs, 
where  the  two  algorithms  for  BEMD  generally  accept  two  types  of  solutions:  (1)  rotating 
signals,  as  intended,  and  (2)  where  the  method  fails  to  extract  the  rotating  components 
and  the  output  signals  are  ones  that  wander  around  zero  in  a  more  complicated  manner 
[11]. 
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Rilling  et  al.  conclude,  from  extensive  simulations,  that  the  outputs  of  the  two 
algorithms  for  BEMD  are  very  similar  when  the  data  clearly  contains  rotating 
components,  but  may  differ  when  they  fail  to  extract  the  rotating  components  [11]. 

2.9.  Comparison  of  EMD  and  Complex  EMD  Extensions 

Yunchao  et  al.  present  the  fact  that  CEMD  is  a  powerful  tool  [15].  The  HHT  is 
the  method  developed  by  Huang  et  al.  [1]  for  analyzing  nonlinear  and  non-stationary 
data,  due  to  the  EMD  technique  not  imposing  any  prior  assumptions  to  the  data  [15]. 

2.9.1.  Realization  of  Complex  EMD 

As  developed  by  Huang  et  al.  [1],  the  original  EMD  method  is  based  on  a 
characteristic  time  scale  defined  by  the  local  extrema.  Yet,  the  original  EMD  method  is 
applicable  only  for  real-valued  time  series  and  it  is  necessary  to  extend  the  application  of 
EMD  into  the  complex  domain.  There  have  been  three  different  methods  proposed  for 
the  realization  of  CEMD  and  while  all  three  have  their  merits,  the  algorithm  proposed  by 
Rilling  et  al.  [11]  is  used  by  Yunchao  et  al.  due  to  the  claim  that  it  is  more  intuitive  [15]. 

2.9.2.  Characteristics  of  Complex  EMD 

From  the  simulations  performed,  the  authors  study  of  the  IMFs  characteristics, 
following  a  method  proposed  by  Zhaohua  et  al.  A  couple  conclusions  that  Yunchao  et  al. 
make  from  these  simulations  were  that:  (1)  CEMD  is  an  effective  dyadic  filter  just  like 
the  original  EMD  method,  (2)  the  power  spectrum  of  the  resulting  complex- valued  IMFs 
are  subject  to  a  normal  distribution,  and  (3)  the  frequency  features  are  the  same  for  the 
real  and  imaginary  part  of  the  same  IMF  [15]. 
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2.9.3.  Numerical  Experiments 

In  addition  to  simulations  on  white  noise,  the  authors  performed  numerical 
experiments  on  a  two-dimensional  (or  velocity  and  pressure)  vector  sensor.  From  these 
numerical  experiments,  Yunchao  et  al.  claim  that  the  decomposed  results  from  the 
CEMD  are  better  than  the  original  EMD  method  [15].  There  are  three  main  differences 
from  the  original  EMD  to  the  CEMD.  First,  the  number  of  IMFs  is  identical  for  the  real 
and  imaginary  parts  with  the  CEMD,  where  the  original  EMD  method  can  have 
superimposed  IMFs,  resulting  in  an  incorrect  number  of  IMFs.  The  second  difference  is 
that,  for  the  CEMD,  the  frequency  characteristic  is  the  same  between  the  same  order 
IMFs  of  the  imaginary  and  real  parts  while  the  original  EMD  method  does  not  possess 
this  characteristic.  This  distinction  between  the  CEMD  and  the  original  EMD  methods  is 
due  to  the  CEMD  observing  the  changes  in  the  two  variables  of  the  signal.  Finally,  it  is 
obvious  from  the  experiments  that  the  characteristics  of  the  IMFs  from  the  analytic  (or 
complex)  signal  are  better  resolved  in  the  frequency  domain  than  the  real  signal  [15]. 

From  the  numerical  experiments  performed  by  Yunchao  et  al.,  it  can  be  concluded 
that,  if  the  signal  is  long  enough  and  the  intricacy  can  be  ignored,  the  analytic  signal  is  a 
good  choice  for  the  application  of  the  CEMD  method;  however,  if  the  signal  is  short  or 
the  system  is  a  real-time  system,  it  is  better  to  choose  the  less  intricate  complex  signal 

[15]. 

Yunchao  et  al.  conclude  that  CEMD  is  consistent  for  the  frequency  characteristics 
of  the  IMFs  for  the  real  and  imaginary  parts.  Additionally,  the  results  decomposed  by  the 
CEMD  method  are  more  legible  than  those  decomposed  by  the  EMD  method,  and  the 
estimations  done  by  the  CEMD  are  more  accurate.  A  final  conclusion  the  authors  make 
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is  that  one  can  choose  the  analytic  signal  or  less  intricate  complex  signal  for  application 
of  the  CEMD  method  in  different  conditions  [15]. 

2.10.  Applications  for  Complex  EMD 

2.10.1.  Multichannel  Information  Fusion 

Mandic  et  al.  wrote  about  information  “fusion”  via  signal  “fission”  in  the 
framework  of  EMD  [2].  The  fission  part  occurs  first  where  the  signal  is  decomposed  into 
its  oscillatory  components,  then  the  fusion  occurs  when  the  IMFs  are  combined  in  an  ad- 
hoc  fashion  to  provide  knowledge  about  the  process  [2].  Mandic  et  al.  claims  that 
extension  of  EMD  into  the  complex  domain  is  especially  important  for  phase-dependent 
processes;  however,  extending  EMD  into  the  complex  domain  is  not  straightforward  and 
depends  heavily  on  the  criterion  for  finding  the  local  extrema  of  the  signal. 

Complex  representation  of  a  signal  can  be  both  intuitive  and  useful  because  the 
amplitude  and  phase  can  be  modeled  simultaneously  [2] .  There  are  many  fields  of  study 
that  use  only  real-valued  data  structures,  but  several  important  signal  processing  areas  use 
complex- valued  data  structures.  EMD  is  a  data  driven  time-frequency  analysis  technique 
that  is  useful  in  the  analysis  of  nonlinear  and  non-stationary  signals. 

One  well-established  information  fusion  model  is  the  waterfall  model.  The 
method  of  EMD  also  performs  both  signal  conditioning  and  feature  extraction,  key 
components  of  the  waterfall  model.  EMD  provides  the  framework  for  unifying 
information  fission  and  fusion;  therefore,  the  aim  of  the  paper  is  to  provide  justification 
for  the  use  of  EMD,  both  real  and  complex,  in  knowledge  extraction  and  information 
fusion  [2]. 
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The  most  intuitive  way  to  extend  EMD  into  the  complex  domain  would  be  to 
apply  EMD  to  the  real  and  imaginary  parts  separately;  however,  in  performing  this,  any 
mutual  information,  such  as  phase,  that  existed  between  the  original  components  is 
ignored  and  lost.  Therefore,  this  paper  examines  the  effectiveness  of  two  of  the  three 
introduced  CEMD  algorithms:  (1)  CEMD  [9]  based  on  the  direct  use  of  the  Hilbert 
transform  properties  and  (2)  RICEMD  [10]  a  generic  expression  of  the  real  EMD. 

Mandic  et  al.  applied  the  CEMD  method  to  data  and  came  up  with  certain 
advantages  and  disadvantages  for  this  method.  The  primary  advantages  that  were  found 
for  this  method  were  that  it  has  a  straightforward  and  intuitive  math  derivation,  acting  as 
a  dyadic  filter  bank.  However,  with  this  method,  the  disadvantages  seem  to  outweigh  the 
advantages.  Not  only  does  this  method  fail  to  reveal  any  synchronized  events  between 
the  data  streams,  but  the  IMFs  are  deprived  of  their  physical  connection  with  the  original 
data  set.  An  ambiguity  exists  at  the  zero  frequency  due  to  the  way  the  math  derivation  is 
formed  and,  finally,  this  method  cannot  be  extended  to  higher  dimensions  due  to  the 
limitation  of  representing  a  signal  by  its  positive  and  negative  frequencies  [2]. 

In  addition  to  applying  the  CEMD  method,  Mandic  et  al.  also  analyzes  the 
RICEMD  method  and  explain  the  advantages  and  disadvantages  inherent  in  the  method. 
The  method  operates  fully  in  the  complex  domain  and  uses  complex  cubic  splines  for 
analysis  of  the  signal  in  the  complex  domain.  Unlike  the  CEMD  method,  the  RICEMD 
method  creates  an  equal  number  of  IMFs  for  the  real  and  imaginary  parts  and  retains  the 
physical  interpretation  of  the  signal.  One  disadvantage  is  the  choice  of  criterion  for 
finding  the  extrema  of  the  complex  signal  is  not  unique,  the  extracted  complex  IMFs  do 
possess  physical  interpretation  [2]. 
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2.10.2.  Single  Vector  Sensor  Application 

Yunchao  et  al.  introduced  CEMD  to  improve  processing  from  a  single  vector 
sensor  of  complex  sound  data  using  HHT  [16].  Yunchao  et  al.  claim  that  CEMD  is  a 
powerful  tool  for  analyzing  complex  data  and  the  results  yielded  in  this  paper  show  that 
CEMD  is  better  in  using  the  information  between  the  correlative  signals.  In  addition,  the 
analytic  signal  is  beneficial  to  direction  estimation  with  different  targets. 

Vector  sensors  are  sensors  that  can  measure  the  pressure  P  and  orthogonal 
components  of  the  particle  velocity,  Vx  and  Vv ,  simultaneously,  and  vector  sensors  can 

also  improve  target  detection  capability.  The  HHT  has  been  utilized  to  identify  the  multi¬ 
targets  using  the  signals  from  the  single  vector  sensor  based  on  the  frequency  feature  of 
the  HHT.  This  use  of  the  HHT  has  led  to  the  development  of  the  Vector  HHT  (VHHT) 
and  the  application  of  the  VHHT  was  analyzed  in  this  paper,  as  well  as  the  improvement 
CEMD  performs  on  the  VHHT  [16]. 

In  two-dimensional  circumstances,  the  vector  sensor  can  simultaneously  measure 
the  pressure  and  orthogonal  components  of  the  velocity,  Vx  and  VY .  From  this  process, 
the  target’s  direction  can  be  obtained  as  follows: 


6  =  arctan 


'pv;' 

KPV'*J 


(36) 


where  the  P  represents  the  pressure  component  of  the  vector  sensor,  and  Vx  and  V 


represent  the  orthogonal  components  of  the  velocity  along  the  x-axis  and  y-axis. 


43 


By  using  the  EMD  to  decompose  the  pressure  and  velocity  components  into  their 
respective  IMFs,  the  real-valued  and  imaginary-valued  IMFs  are  obtained.  The  analytic 
signals  of  the  resulting  IMFs  are  represented  below: 


ZPj  (t)  =  APj  (t)  exp  (  j(pPj  (0)  (37) 

ZyxJ  (0  =  Ay  (0  exp  (  (tj)  (38) 

? 

and 

zVyj  (0  =  Ayj  (0  exP  ( JAj  (0) 

The  instantaneous  sound  energy  flows  corresponding  to 
represented  by  the  multiplication  of  the  complex  conjugate  of 
analytic  signal  and  the  pressure  analytic  signal, 

SPVJ(t)  =  (zpj(tj)(zvJ(t))  (40) 

and 

S„,ft)=(ZP1(»)(zr;(»)  (41) 


(39) 

the  x-axis  and  y-axis  are 
the  velocity  components 
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Finally,  the  instantaneous  azimuth  of  the  acoustic  signal  is  represented  by  the  following 


equation: 


Oj  it)  =  arctan 


"Rej 

Re 

{Spvjto}  ^ 

(42) 


CEMD  is  an  extension  of  EMD  in  the  complex  plane.  The  method  employed  in 
this  paper  was  the  method  proposed  by  Rilling  et  al.  [11].  In  this  article,  a  number  of 
characteristics  of  the  CEMD  method  are  discovered.  Just  a  few  of  those  characteristics 
include:  (1)  that  CEMD  is  shown  to  act  as  a  dyadic  filter,  (2)  that  the  period  of  the  IMF 
increases  when  the  order  increases  and  the  center  frequency  decreases,  (3)  the  frequency 
feature  of  the  real  and  imaginary  parts  of  the  same  IMF  are  the  same,  and  (4)  that  the 
CEMD  algorithm  is  adaptive  like  the  original  EMD  algorithm  [16]. 

Through  the  experiments  performed,  Yunchao  et  al.  concluded  that  CEMD  takes 
full  advantage  of  the  information  between  the  relevant  signals  in  the  acoustic  vector 
sensors.  There  are  three  separate  aspects  of  the  findings  that  show  how  well  the  CEMD 
works  in  the  field  of  vector  sensors.  First,  the  order  and  frequency  feature  of  the  real  and 
imaginary  parts  of  the  IMFs  are  identical.  Second,  the  CEMD  is  better  than  EMD  in 
noise  suppression  and  reducing  the  mode-mixing.  And  last,  that  the  analytic  signal  is 
more  suitable  for  high  signal-to-noise  (SNR)  and  the  simple  complex  signal  is  better  in 


low  SNR  [16]. 


2.10.3.  Multiscale  Image  Fusion  Application 

Looney  and  Mandic  propose  a  solution  to  the  problem  of  uniqueness  when 
performing  fusion  of  data  from  multiple  and  heterogeneous  sources  [17].  The  proposed 
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solution  relies  heavily  on  using  complex  extensions  of  the  data-driven  technique  of  EMD, 
a  new  analysis  technique  proposed  by  Huang  et  al.  [1]. 

In  the  data  and  information  fusion  community,  there  is  a  significant  challenge 
when  different  focal  points  are  observed  [17].  The  technique  of  EMD  has  been  proposed 
for  data  fusion,  where  only  the  “relevant”  IMFs  are  recombined  into  a  restored  signal. 
Due  to  its  adaptivity,  it  is  natural  to  consider  the  use  of  the  EMD  method  for  the  problem 
of  heterogeneous  image  fusion;  however,  there  is  still  a  problem  of  uniqueness  when 
using  the  EMD  method.  Therefore,  extensions  into  the  complex  domain  of  the  EMD 
method  have  been  proposed  to  help  with  the  uniqueness  of  the  resulting  IMFs  [17]. 
While  there  have  been  three  complex  extension  of  EMD  recently  proposed,  the  method 
employed  in  this  paper  BEMD  [11]. 

The  local  and  data-driven  nature  of  EMD  leads  to  two  problems.  First,  that  the 
uniqueness  of  decomposition-signals  gives  different  IMFs,  and  second,  that  mode-mixing 
of  the  IMFs  occurs  [17].  The  problem  of  uniqueness  can  be  addressed  by  stopping  the 
decomposition  once  a  specific  number  of  IMFs  has  been  obtained.  If  the  number  of  IMFs 
from  each  source  is  equal  in  number,  then  the  problem  of  mode-mixing  is  also  fixed.  So, 
in  order  to  fix  both  problems,  Looney  and  Mandic  propose  to  apply  the  BEMD  method  to 
decompose  heterogeneous  complex  data  simultaneously,  rather  than  decomposing  one 
part  of  the  signal  at  a  time  [17]. 

To  show  the  effectiveness  of  using  EMD,  and  specifically  BEMD,  simulations 
were  performed  by  Looney  and  Mandic  on  generated  complex  data  and  real-world  fusion 
data.  An  automatic  fusion  algorithm  is  also  proposed  that  is  based  on  the  BEMD 
algorithm  [17].  The  robustness  of  the  analysis  of  the  generated  data  guarantees  a 
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meaningful  comparison  between  the  scales  and  forms  the  basis  for  the  proposed  image 
fusion  algorithm  and  will  be  used  for  the  real-world  data  analysis.  By  applying  the 
BEMD  algorithm  to  real-world  image  fusion  data,  it  is  illustrated  that  the  problems  of 
mode-mixing  and  uniqueness  can  be  easily  overcome  [17]. 

Looney  and  Mandic  conclude  that  the  potential  for  BEMD  for  information  fusion 
is  verified,  as  well  as  a  set  of  common  frequency  scales  can  be  determined  by 
simultaneously  decomposing  sources  using  the  BEMD  method  [17].  In  addition,  the 
application  of  BEMD  enables  the  proposed  approach  to  overcome  the  uniqueness  and 
mode-mixing  problems.  For  future  work,  Looney  and  Mandic  propose  that  higher 
dimensional  extensions  should  be  developed  in  order  to  enable  the  fusion  of  more  than 
two  images  [17]. 

2.11.  Multivariate  EMD  Application 

Mutlu  and  Aviyente  present  the  importance  of  quantifying  the  phase  synchrony 
between  signals  is  stated  for  different  applications  [12].  However,  current  techniques 
used  to  measure  and  quantify  the  phase  synchrony  suffer  from  constraints  inherent  to  the 
wavelet  transform  and  Hilbert  transform.  Therefore,  in  order  to  address  such  constraints 
on  the  analysis  of  the  signals,  a  recently  introduced  multivariate  empirical  mode 
decomposition  (MEMD)  in  order  to  assist  in  the  quantification  of  multivariate  phase 
synchrony. 

Mutlu  and  Aviyente  propose  to  use  MEMD  for  quantifying  the  phase  synchrony 
between  multiple  time  series  [12].  The  original  EMD  method  acts  as  a  dyadic  filter,  thus 
a  “pre-filtering  tool”  for  the  Hilbert  transform-based  phase  synchrony  analysis.  The  goal 
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of  Mutlu  and  Aviyente’s  research  is  to  extend  the  measures  of  correlation  for  multiple 
variables  from  statistics  for  quantifying  multivariate  synchronization.  Mutlu  and 
Aviyente  claim:  (1)  the  MEMD  might  be  used  to  define  pairwise  synchrony  between 
multiple  time  series  across  the  same  frequency,  and  (2)  MEMD  can  extend  the  notion  of 
bivariate  synchrony  to  multivariate  synchronization  [12]. 

From  their  research,  a  new  approach  for  quantifying  multivariate  phase 
synchronization  within  a  group  of  oscillators.  The  new  approach  is  based  on  the 
application  of  MEMD  for  extracting  time  and  frequency  dependent  phase  information 
[12].  MEMD  method  results  in  two  improvements  discovered  from  their  research.  First, 
MEMD  is  data-driven  and  eliminates  the  need  for  chosen  bandpass  filters  and  second,  the 
MEMD  extends  the  current  state  of  the  art-phase  synchrony  analysis  from  quantified 
bivariate  relationships  to  the  multivariate  case  [12].  Mutlu  and  Aviyente  propose  for 
future  work  that  focuses  on  the  extension  of  the  methods  proposed  using  different 
multivariate  analysis  technique. 

2.12.  Conclusions 

Rilling  et  al.  conclude  that  the  EMD  technique  is  a  promising  tool  but  it  needs  to 
be  better  understood.  They  call  for  further  studies  devoted  to  a  theoretical  approach  and 
closed-form  solution,  due  to  the  EMD  definition  by  an  algorithm  rather  than  a  closed- 
form  solution.  Additionally,  Kim  and  Oh  conclude  that  the  IMFs  that  are  decomposed 
from  the  EMD  technique  “provide  a  multi-resolution  tool  and  spectral  analysis  given 
local  information  with  time- varying  amplitude  and  phase”  [13]. 
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The  CEMD  method  takes  full  advantage  of  the  information  between  the  relevant 
signals  in  the  acoustic  vector  sensors  [16].  In  addition,  the  results  decomposed  by  the 
CEMD  method  are  more  legible  than  those  decomposed  by  the  EMD  method,  and  the 
estimations  done  by  the  CEMD  are  always  more  accurate  [15]. 

Concerning  the  three  CEMD  methods,  conclusions  have  been  made  about  the 
possibilities  inherent  to  each  one  and  the  disadvantages  for  each  method.  The  first 
complex  method,  ambiguously  called  CEMD,  presented  conclusions  that  the  CEMD 
method  can  be  achieved  based  on  some  inherent  properties  of  the  complex  signals  [9]. 
For  the  RICEMD  method,  Umair  Bin  Altaf  et  al.  conclude  that  the  analysis  of  real-world 
complex-valued  data  shows  that  the  proposed  method  provides  new  insights  into  time- 
frequency  analysis  of  nonlinear  and  non- stationary  signals,  which  was  not  possible  before 
[10].  The  authors  of  the  BEMD  method  conclude  that  the  outputs  of  the  two  algorithms 
for  BEMD  are  very  similar  when  the  data  clearly  contains  rotating  components,  but  may 
differ  when  they  fail  to  extract  the  rotating  components  [11]. 

Through  all  the  papers  written  for  the  EMD  technique,  each  application  of  the 
method  shows  that  the  method  proves  worthwhile  when  applied  to  the  data  sets.  There 
are  some  restraints  still  inherent  to  the  original  EMD  method;  various  authors  then  extend 
the  method  into  the  complex  domain,  which  helps  with  some  of  the  issues  found  in  the 
real  domain.  Overall,  the  investigation  into  the  HHT  developed  by  Huang  et  al.  [1],  and 
more  specifically  the  EMD  technique,  displays  information  not  accessible  by  the  FFT  and 
that  new  information  can  help  in  the  analysis  of  many  natural  phenomena  that  are 
nonlinear  and  non- stationary  in  nature. 
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III.  Algorithm  Design  and  Implementation 


3.1.  Overview 

This  chapter  presents  the  methodology  used  in  achieving  the  objective  of  the 
research.  The  objective  of  this  research  was  two-fold.  First,  the  assessment  of  using  the 
CEMD  technique  inherent  to  the  HHT  as  a  signal  processing  analysis  tool  was 
investigated.  Secondly,  the  results  from  the  application  of  the  Hilbert  transform  and  the 
FFT  on  the  decomposed  data  sets  using  the  CEMD  technique  were  compared.  The 
comparisons  of  the  resulting  graphs  from  the  application  of  the  two  separate 
mathematical  transforms  were  sought  to  deliver  insight  into  determining  which  of  the 
transforms  provided  enhanced  fidelity  of  the  real-world  data  set.  Additionally,  the 
analysis  of  the  HHT  was  compared  with  the  analysis  provided  from  the  FFT  to  determine 
how  the  application  of  the  HHT  affects  a  RCS  data  set. 

The  approach  used  to  satisfy  the  two  requirements  stated  above  concerning  the 
problem  statement,  was  conducted  primarily  using  the  MATLAB®  computer  program, 
where  the  Signal  Processing  and  Spline  Toolboxes  were  employed.  In  addition  to  using 
MATLAB®  and  the  two  above  mentioned  toolboxes,  code  that  was  written  by  Dr. 
Flandrin  and  available  on  his  website  was  implemented,  as  detailed  in  BEMD  [11]. 

The  method  for  collection  of  the  data  will  be  discussed,  as  well  as  the  components 
of  the  designed  algorithm  will  be  presented  and  explained.  Finally,  a  detailed 
explanation  of  the  implementation  of  the  algorithm  used  will  be  discussed,  with  detailed 
explanations  of  the  various  components  used  in  completion  of  the  methodology. 
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3.2.  Data  Collection 


There  were  two  separate  data  collections  employed  for  the  completion  of  this 
evaluation.  The  first  data  set  was  provided  by  the  sponsor,  also  referred  to  as  the  “real- 
world  data  set”  during  the  methodology  portion  of  this  thesis.  This  set  of  data  was 
obtained  through  determining  the  complex  RCS  data  from  a  rotating  object. 
Measurement  of  a  target's  RCS  is  performed  at  a  radar  reflectivity  range  or  scattering 
range.  One  type  of  RCS  range  is  an  outdoor  range,  where  the  target  is  positioned  on  a 
specially  shaped  low  RCS  pylon  some  distance  down-range  from  the  transmitters.  Such  a 
range  eliminates  the  need  for  placing  radar  absorbers  behind  the  target;  however,  multi- 
path  interactions  with  the  ground  must  be  mitigated.  There  are  instances,  as  with  the 
real-world  data  set,  in  which  an  object  exists  in  the  atmosphere  and  pulses  are  sent  from 
the  radar  to  the  object,  recording  the  radar  return  from  the  object.  The  real-world  data  set 
collected  consisted  of  at  least  two  revolutions  of  the  object,  while  the  second  set  of  data 
consisted  of  only  one  extremely-finely  sampled  revolution  of  the  object.  For  the  real- 
world  object,  the  approximate  shape  was  a  cylinder  with  one  end  enclosed,  with  the  other 
end  open,  or  a  “cavity”  cylinder. 

The  second  data  collected  consisted  of  using  the  RCS  range  owned  by  AFIT 
called  an  anechoic  chamber.  In  such  a  room,  the  target  is  placed  on  a  rotating  pillar  in  the 
center,  and  the  entire  background  is  covered  with  radar  absorbing  material.  These 
absorbers  prevent  corruption  of  the  measurement  due  to  reflections.  A  compact  range  is 
an  anechoic  chamber  with  a  reflector  to  simulate  far-field  conditions.  There  were  a 
variety  of  shapes  used  over  the  course  of  the  data  collection,  but  upon  further  analysis  of 
the  real-world  object,  four  specific  sets  of  data  were  chosen  for  comparison  with  the 
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decomposed  real-world  data  set.  These  four  data  sets  were:  (1)  a  cylinder  with  two  end 


caps,  refer  to  Figure  1;  (2)  a  cylinder  with  one  end  cap  and  one  open  “cavity”,  refer  to 
Figure  2;  (3)  a  cone-sphere,  or  an  object  with  a  roundly  tapered  end,  refer  to  Figure  3;  and 
(4)  a  dihedral  corner  reflector  configuration,  refer  to  Figure  4. 


Figure  1:  Cylinder  with  two  end  caps,  on  Styrofoam  pylon 


Figure  2:  Cylinder  with  one  end  cap  and  one  open  end,  “cavity”,  on  Styrofoam  pylon 


Figure  3:  Cone-Sphere,  on  Styrofoam  pylon 
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Figure  4:  Dihedral  Spherical  Corner  Reflector,  on  Styrofoam  pylon 

The  four  data  sets  detailed  above  were  considered  the  simulated  (or  theoretical) 
data  to  help  assess  whether  the  CEMD  technique  worked  properly  for  a  signal  composed 
of  various  canonical  objects,  such  as  the  dynamic  object  represented  by  the  real-world 
data  set.  To  assess  the  usefulness  and  workability  of  the  CEMD  algorithm,  the  real-world 
data  set  was  decomposed  and  compared  with  the  RCS  values  of  the  simulated  objects  to 
determine  whether  those  components  were  possibly  detected  in  the  decomposed  real- 
world  signal. 

To  work  with  a  complex- valued  data  set  in  the  same  manner  as  the  real-world 
data  set,  the  simulated  data  collected  in  the  anechoic  chamber  had  to  be  calibrated.  This 
calibration  was  accomplished  using  a  MATLAB®  GUI  called  ALPINE©,  version  3.1.1, 
designed  by  AFIT  professor  Dr.  Peter  Collins  [18].  Through  the  use  of  two  specific 
modules  of  the  ALPINE©  GUI,  the  calibrateRCS  and  the  plotGlobalRCS  modules,  the 
data  input  and  calibration  was  user-friendly  and  accomplished  quickly.  Access  to  this 
GUI  saved  hours  of  individual  code  creation  and  analysis.  The  calibration  results  were 
stored  in  a  MATLAB®  MAT-file,  enabling  the  flexibility  to  load  the  data  into  any 
desired  program. 
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Once  the  data  was  calibrated  and  placed  in  a  useable  format,  the  data  was  able  to 
be  analyzed  using  the  CEMD  algorithm  developed  by  Dr.  Flandrin  [11],  followed  by  the 
application  of  both  the  FFT  and  the  Hilbert  transform  to  the  IMFs  decomposed  using  the 
CEMD  algorithm. 

3.3.  Algorithm  Design  and  Components 

Traditionally,  signal  processing  data,  such  as  the  RCS  data  collected  for  this 
thesis,  was  analyzed  using  Fourier  spectrum  analysis  [1].  However,  the  technique  of  the 
EMD  method  and  use  of  the  complete  HHT  have  been  receiving  more  attention  for  use  in 
the  analysis  of  signal  processing  data.  More  specifically,  the  CEMD  algorithm  was 
employed  for  the  decomposition  of  the  signals  analyzed  in  this  thesis.  In  this  section,  the 
original  EMD  algorithm  was  explored,  as  well  as  the  complex  extension  of  the  EMD 
algorithm.  Once  these  two  items  have  been  discussed,  explanation  of  the  components  in 
the  algorithm  employed  for  the  thesis  will  be  detailed. 

3.3.1.  The  Empirical  Mode  Decomposition 

EMD  is  defined  as  an  exploratory  analysis  technique.  EMD  is  an  adaptive 
technique  used  to  decompose  a  given  signal  into  its  oscillatory  modes  [2].  This 
decomposition  is  accomplished  through  a  process  referred  to  as  the  sifting  algorithm. 
The  sifting  algorithm,  which  defines  the  EMD  process,  was  detailed  previously  in  Table 
1  [1]- 

Once  the  sifting  algorithm  was  applied  to  the  given  signal,  the  resulting 
oscillatory  components  are  called  IMFs  [2].  These  IMFs  represent  the  oscillatory  nature 
embedded  in  the  data.  In  addition,  these  IMFs  also  represent  the  basis  functions  of  the 
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signal  and  are  derived  from  the  data,  rather  than  a  pre-defined  set  of  basis  functions  that 
exist  in  such  transforms  as  Fourier  and  wavelets  [1], 

The  algorithm  detailed  previously  in  Table  1  was  utilized  to  perform  the  EMD 
decomposition  as  follows: 

x(k)  =  ^=1ci(k)  +  r(lc)  (43) 

in  which  the  IMFs  are  denoted  by  ct(k)  and  the  residue  is  represented  by  r(k ) . 

Once  the  sifting  algorithm  was  applied  to  the  given  signal,  the  IMFs  are  in  a  form 
that  can  be  linearly  transformed.  More  specifically,  the  FFT  and  Hilbert  transform  were 
applied  to  the  IMFs  of  the  decomposed  data.  Through  application  of  the  Hilbert 
transform,  the  given  real-valued  signal  was  transformed  into  a  complex-valued  signal, 
where  the  analytic  representation  of  (43)  is  given  by: 


X(t)  =  ^a,(t)em,>  (44) 

In  (44),  the  residue  r(t)  was  omitted  due  to  its  lack  of  oscillatory  behavior  [2]. 
The  analytic  signal  was  created  using  the  IMF  to  represent  the  real  part  and  employing 
the  Hilbert  transform  to  represent  the  imaginary  part,  such  that  x+  /T I  (a  )  becomes  the 
new  representation  of  the  signal  [2].  The  aim  for  creating  this  analytic,  time -dependent 
signal  was  to  extract  the  time-dependent  amplitude  a;  (/)  and  the  phase  function,  0j  (/) , 
components  more  easily  [2] .  Finally,  the  instantaneous  frequency,  denoted  by, 


co(t) 


dO(t) 
dt  , 


(45) 


can  also  be  extracted,  used  to  create  a  Time-Frequency-Amplitude  representation  of  the 
signal  called  the  Hilbert  spectrum  [2]. 
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3.3.2.  Complex  Extension  of  Empirical  Mode  Decomposition 
The  original  formulation  of  the  EMD  algorithm  by  Huang  et  al.  [1]  restricts  its 
application  to  real-valued  signals.  However,  many  fields  of  study  use  complex-valued 
data  sets,  whereby  ignoring  the  imaginary  components  can  cause  important  information 
to  be  lost,  as  well  as  the  analysis  of  the  signal  incompletely  presented.  Since  introduction 
of  the  EMD  algorithm  by  Huang  et  al.  in  1998,  there  have  been  three  complex  extensions 
proposed  [9,10,11].  In  Chapter  II  section  2.7,  more  detailed  descriptions  of  the  CEMD 
[9]  and  the  RICEMD  [10]  methods  are  located.  The  third  complex  extension.  BEMD 
[11],  was  revisited  due  to  its  involvement  in  the  research  and  the  decomposition  of  given 
complex  RCS  data  sets. 

3.3.2. 1  Bivariate  Empirical  Mode  Decomposition 
The  method  proposed  by  Rilling  et  al.  [11]  was  followed  to  study  the 
characteristics  of  the  IMFs  of  the  given  data  sets.  The  other  two  complex  extensions 
discussed  in  Chapter  II  employed  the  original  EMD  algorithm  and  decomposed  complex¬ 
valued  signals  similar  to  real- valued  signals;  however,  the  BEMD  algorithm  “adapts  the 
rational  underlying  the  EMD  to  the  bivariate  [or  complex-valued]  framework”  [11].  Dr. 
Flandrin’s  BEMD  paper  provided  free  MATLAB®  code  [11]  for  implementation  of  the 
CEMD  algorithm,  Hilbert  transform,  and  Fast-Fourier  transform. 

The  original  EMD  algorithm  is  based  on  the  natural  oscillation  related  to  the 
signal  extrema.  Nevertheless,  for  a  complex-valued  signal,  defining  extrema  in  the  same 
manner  as  a  real- valued  signal  is  more  confusing  and  unclear.  Therefore,  Rilling  et  al. 
proposed  the  notion  of  “rotation”  and  is  considered  a  three-dimensional  interpretation  of 
the  real-valued  notion  of  oscillation  [11].  The  underlying  idea  of  the  complex  extension 
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to  EMD  proposed  by  Rilling  et  al.  is  that  the  “bivariate  signal  =  fast  rotations 
superimposed  on  slower  rotation”  [11],  whereas  the  original  EMD  algorithm  is  based 
upon  the  idea  that  the  “signal  =  fast  oscillations  superimposed  on  slower  oscillations” 
[11].  As  with  the  original  EMD  algorithm,  the  BEMD  algorithm  uses  the  idea  of  an 
“envelope,”  where  the  envelope  of  the  BEMD  algorithm  is  now  a  three-dimensional  tube 
enclosing  the  signal  [11]  (refer  to  Figure  5). 


Figure  5:  The  signal  enclosed  in  its  3D  envelope.  The  black  thick  lines  stand  for  the  envelope  curves 

that  are  used  to  derive  the  mean. 

Using  this  idea  of  how  the  envelope  is  created,  the  slowest  rotating  component  is 
defined  as  the  center  of  the  tube.  There  are  two  ways  of  defining  the  mean  value  of  the 
tube:  (1)  define  the  mean  as  the  barycenter  (or  center  of  mass)  of  the  four  points 
considering  each  to  have  unit  mass  (refer  to  Figure  6),  or  (2)  define  the  mean  as  the 
intersection  of  two  straight  lines,  one  being  halfway  between  the  two  horizontal  tangents 
and  the  other  one  halfway  between  the  two  vertical  tangents  (refer  to  Figure  7)  [11]. 
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Figure  6:  Illustration  of  the  first  definition  of  the  mean  of  the  complex-valued  signal. 


Figure  7:  Illustration  of  the  second  definition  of  the  mean  of  the  complex- valued  signal. 

Due  to  its  natural  robustness  to  sampling  errors,  scheme  two  for  defining  the 
mean  value  of  the  tube  is  preferred  in  practice.  These  sampling  errors  should  be  taken 
seriously  since  the  original  EMD  is  sensitive  to  sampling,  leading  to  the  idea  that  the 
BEMD  algorithm  will  also  be  sampling  sensitive  [11]. 

The  desired  interpolation  is  defined  as:  “a  smooth  interpolation,  with  as  few 
‘spurious  bumps’  as  possible”  [10].  To  satisfy  this  definition,  a  cubic  spline  was 
employed.  The  interpolation  of  the  BEMD  algorithm  is  performed  in  a  similar  manner  to 
the  EMD  algorithm.  Wolfram -MathWorldl  defines  a  cubic  spline  as: 
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“A  spline  constructed  of  piecewise  third-order  polynomials  which  pass  through  a 
set  of  m  control  points.  The  second  derivative  of  each  polynomial  is  commonly 
set  to  zero  at  the  endpoints,  since  this  provides  a  boundary  condition  that 
completes  the  system  of  m  —  2  equations.  This  produces  a  so-called  ‘natural’ 
cubic  spline  and  leads  to  a  simple  tri-diagonal  system  which  can  be  solved  easily 
to  give  the  coefficients  of  the  polynomials”  [19]. 

The  proposed  complex  extension  of  the  BEMD  algorithm  led  to  the  complex 
EMD  algorithm  detailed  previously  in  Tables  4  and  5.  The  proposed  algorithms  used  the 
same  algorithm  as  the  EMD,  where  the  only  difference  was  the  new  sifting  operators, 
SBl  and  SB2 ,  representing  the  fast  and  slow  oscillations. 

The  reformulation  of  the  second  complex  algorithm  allowed  the  sifting  operator 
to  be  represented  as  a  univariate  (or  real-valued)  EMD  sifting  operator,  shown  below  in 
Table  6.  This  reformulation  allowed  the  behavior  of  the  algorithm  to  be  studied  in  a 
similar  method  as  the  original  EMD  algorithm  [11]. 


Table  6:  Algorithm  2  Reformulation  for  EMD  Bivariate  Extension 


2  kn 


1.  Given  an  angular  direction:  pk  = -  for  1  <  k  <  N  complete  steps  2  through  3 


N 


2.  Project  the  complex-valued  signal  x(t)  on  direction  pk:dp/_  =R e,[eJPkx(tj\ 


3.  Compute  the  partial  estimate  in  direction  pk  :  spj  =e  ;aP[<9a  ](f) 


Pk  ' 


4.  Subtract  the  mean  to  obtain  SB2(t)  -  —  V  v  ( t ) 

v  '  yy  Z— Ik  Pk  v  ' 


3. 3. 2. 2.  Bivariate  Intrinsic  Mode  Functions  (IMFs) 

The  BEMD  algorithm  proposed  by  Rilling  et  al.  was  designed  so  “that  signals 
rotating  around  zero  are  admissible  outputs”  [11].  Rilling  et  al.  clarified  this  vague 
notion  of  “rotating  around  zero”  by  asking:  “what  signals  the  algorithms  actually  consider 
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admissible  outputs”  [11].  In  other  words,  x(t )  is  considered  a  “fixed  point  of  the  sifting 
operator”  [1 1],  or 


SB1x(0  «  x(t ) 


SB2x(t)  ~  x(t)  ■ 


(46,  47) 


To  explain  the  outputs  that  resulted  from  the  BEMD  algorithm,  Rilling  et  al. 
explained  various  simulations.  Both  algorithms  presented  by  Rilling  et  al.  state  two  types 
of  solutions  are  generally  accepted.  The  first  solution  corresponded  to  the  rotating 
signals,  as  expected,  and  the  second  occurred  in  cases  where  the  method  fails  to  extract 
the  rotating  components.  In  the  latter  case,  the  outputs  “wander  around  zero”  was  more 
complicated  than  the  first  case  [11].  The  second  types  of  solutions  were  encountered 
when  the  signal  was  not  a  clearly  rotating  signal,  whereas  solutions  of  the  first  type  were 
signals  with  unchanging  local  rotation  [11]. 

Rotating  “around  zero,”  as  stated  by  Rilling  et  al.,  was  clarified  by  considering  the 
simple  case  where  the  signal  performs  one  rotation  around  zero  per  period,  or 


i//(t  +  T)  =  i//(t)  +  2n ,  (48) 

where  T  is  the  period  of  the  signal.  The  clarification  allows  the  understanding  that  only 
one  maximum  value  exists  per  period.  Thus,  all  of  the  envelope  curves  are  constants 
with  respect  to  time,  allowing  the  mean  to  be  derived  analytically.  Consequently,  the 

envelope  curve  associated  with  the  angular  direction,  pk  =  ,  is  equal  to  the 

maximum  signal  value  in  that  specific  direction,  where  the  phase  of  the  derivative  is 
equal  to  the  definition  in  (48). 
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For  the  two  separate  algorithms  existing  in  the  BEMD  method,  Rilling  et  al. 
showed  the  means  are  equivalent  to  each  other  [11].  The  mean  for  the  first  algorithm 
presented  in  Table  4  is  defined  as: 


m 


\t)  =  —  V  x 


( 

-1 

r  2  7tk 

V 

+— 

V 

l  N 

2  )j 

with  the  limit  of  the  algorithm  resulting  in: 


(49) 


( p))dp 


(50) 


representing  the  mean  of  the  signal  over  a  period  with  a  weight  function  of 


dy/ 

dt 


>0.  The 


weighting  conveys  that  the  sampling  is  denser  where  the  curvature  is  larger. 

In  a  similar  manner,  the  mean  for  the  second  algorithm  or  limit  notation,  omitting 
the  summation  notation,  is: 


mB2(t)  =  -  T em,)K {e~m,)x{t)}^-dt 
n Jo  dt 

=  mm(t)  +  —  fVy(V(0^  dt 

2  7t]o 


dt 


This  equation  leads  to  mB2(t )  =  mB\t)  due  to: 


B 1  . 


(51) 


?ejl//(>)x*  (t}  — dt  = 


dy/ 

dt 


e~m,)x 


(,)T  +  r^e^'d, 

JO  Jo  Af 


0  dt 

=  ^\T  r(t)em,)dt  =  l\l—dt=  0  . 

2  Jo  ?  Jo  rtf 


2  Jo  dt 


(52) 


The  mean  is  equivalent  for  both  BEMD  algorithms  and  is  also  a  fixed  point  of 
both  sifting  operators  if  and  only  if  the  integral  in  (50)  is  close  to  zero.  More  generally, 
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the  outputs  of  the  two  complex  algorithms  are  similar  when  there  are  rotating 
components  inherent  to  the  complex  data,  but  may  differ  when  the  algorithms  fail  to 
extract  the  rotations. 

3.3.3.  The  Hilbert  Transform 

A  real  function,  fit) ,  and  its  Hilbert  transform,  fit) ,  are  related  where,  together, 
they  create  an  analytic  signal.  This  analytic  signal  can  be  represented  as  amplitude  and 
phase,  with  the  derivative  of  the  phase  called  the  instantaneous  frequency  [20].  By  taking 
the  Fourier  transform  of  a  strong  analytic  signal  as  described  above,  the  “negative” 
frequencies  are  discarded  and  the  spectrum  becomes  one-sided  in  the  frequency  domain. 

The  Hilbert  transform  in  the  time  domain  is  a  convolution  between  the  Hilbert 

transformer,  ,  and  the  function,  f(t) .  The  Hilbert  transform,  H\  f(t)\ ,  of  f(t) 
defined  as: 

H[f(t)]  =  fit )*  —  =  -f  —  dr  =  -  f  .  (53) 

7tt  7T:^°t-T  U T 

This  convolution  represents  the  response  to  fit)  of  a  linear  time  invariant  filter 
having  the  impulse  response,  [21].  The  Hilbert  transform  was  expressed  as: 

-  1  r°°  f(r) 

H[fit)]  =  fit)  =  -P\  — dr .  (54) 

The  integral  in  the  Hilbert  transform  description  in  (54)  is,  by  definition, 
improper.  The  integrand  contains  a  singularity,  but  has  infinite  limits  of  integration  [21]. 
The  P  in  front  of  the  integral  represents  the  Cauchy  principal  value,  defined  by  the 
following  equation: 
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=  —  lim  f  r^dT+ry^dT 

ne^  \h-yet-T  it+e  r_r 


(55) 


The  Cauchy  principal  value  was  obtained  by  considering  a  finite  range  of 
integrate,  symmetric  about  the  singularity,  if  and  when  it  exists  [21]. 

3.3.3. 1  Mathematical  Motivations 

The  signal  has  a  Fourier  transform  of: 


-j  sgn(<y)  = 


-j,ifa»  0 
0  ,ifco  =  0 
j,ifox  0 


(56) 


If  the  signal  f(t )  has  the  Fourier  transform  of  F(oj)  ,  then  the  Hilbert  transform,  fit) , 


has  the  Fourier  transform  of  [21] 


F(co)  =  -j  sgn (co)F(co) .  (57) 

The  Hilbert  transform  is  more  easily  understood  in  the  frequency  domain.  The 
magnitude  of  F(a>)  does  not  change;  however,  the  phase  of  F(co)  does  change.  The 
positive  frequencies  of  the  Fourier  transform  values  are  multiplied  by  -j  (or  a  phase 

change  of  ),  while  the  negative  frequencies  of  the  Fourier  transform  values  are 

multiplied  by  j  (or  a  phase  change  of  ).  Therefore,  given  the  signal  spectrum 

F(o))  -a  +  jb ,  the  Hilbert  transform  is  F(co)  =  b-  ja  for  co>  0  and  F(co)--b  +  ja  for 
0)<  0 .  The  Hilbert  transform  exchanges  the  real  and  imaginary  parts  of  the  signal 
Fico)  =  a  +  jb ,  while  changing  one  of  the  signs  in  the  signal  [21]. 


63 


A  table  of  common  Hilbert  transform  pairs  is  detailed  below  in  Table  7  [21]. 


Table  7:  Hilbert  Transform  Pairs 


Signal,  fit) 

Hilbert  transform,  /  (t) 

algl(t)  +  a2g2(t);a1,a2  eC 

«iii(0+«2i2(0 

h(t-t0 ) 

o 

i 

h(at)',a  ^  0 

sgn  {a)hiat) 

d  ,  .  x 

—hit) 

dt 

d  ; ,  . 

—  hit) 

dt 

S(t) 

1 

Kt 

eJl 

~je" 

e-jt 

je  " 

COS  it) 

sin(/) 

rectit) 

1,  2t  +  l 

—  In - 

n  2t-l 

sin  cit) 

Ttt  .  ^  \  .  ( Ttt 

—  sinc“(r)  =  sin  — 

2  V  2 

\ .  fO 

sine  — 

)  v2  J 

1 

1  +  r 

t 

1  +  t2 

The  Fourier  transform  is  important  for  signal  processing.  When  given  a  real 


function  fit) ,  only  the  positive  frequency  axis  is  of  interest  because  spectrum  is 
symmetric  about  zero.  The  negative  frequency  axis  is  not  needed  and  the  Hilbert 
transform  is  employed  to  remove  the  negative  frequency  axis  [20]. 

The  Hilbert  transform  of  a  strong  analytic  signal  is: 


=  H[f(t)  +  jf(t )]  =  /(*)-  jfit)  =  -jz.(t) .  (58) 

The  Hilbert  transform  can  be  used  to  create  an  analytic  signal  from  a  real  signal 

[20], 
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Therefore,  it  is  possible  to  study  the  signal  as  a  rotating  vector  with  an 
instantaneous  phase  pit )  and  instantaneous  amplitude  A(l)  in  the  time  domain: 


z(t)  =  f(t)  +  jf(t)  =  A(t)ei«,).  (59) 

The  final  portion  of  (59)  represents  the  signal  in  polar  notation,  where 

=  +  (60) 

and 


6(t )  =  arctan 


Finally,  the  notion  of  instantaneous  frequency  is  introduced  as: 


(61) 


co(t) 


dO(t) 

dt 


3. 3. 3. 2  MATLAB®  Hilbert  Transform  Function 


(62) 


MATLAB®  has  a  function  that  will  perform  the  Hilbert  transform  as  defined  in 
this  section.  Under  the  help  section  for  the  hilberto  function  in  MATLAB®,  the 
explanation  of  the  input  values  and  the  outputs  are  described  in  great  detail.  The 
hilberto  function  returns  a  complex  sequence,  or  the  analytic  signal,  from  a  real¬ 
valued  sequence.  The  real-valued  part  of  the  signal  contains  the  original  signal  of  the 
input  data  and  the  imaginary-valued  part  of  the  signal  contains  the  Hilbert  transform  of 
the  input  data  [22].  The  imaginary-valued  portion  represents  a  90  degree  phase  shift  of 
the  original  real-valued  data  sequence.  In  addition,  “the  Hilbert  transformed  series  has 
the  same  amplitude  and  frequency  content  as  the  original  real  data  and  includes  phase 
information  that  depends  on  the  phase  of  the  original  data”  [22]. 
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MATLAB®  details  the  algorithm  inherent  to  the  MATLAB®  hilbert ( ) 
function,  as  the  following: 

“The  analytic  signal  for  a  sequence  x  has  a  one-sided  Fourier  transform,  that  is, 
negative  frequencies  are  0.  To  approximate  the  analytic  signal,  hilbert 
calculates  the  FFT  of  the  input  sequence,  replaces  those  FFT  coefficients  that 
correspond  to  negative  frequencies  with  zeros,  and  calculates  the  inverse  FFT  of 
the  result”  [22]. 

In  greater  detail,  the  hilbert  ()  function  employs  the  following  four-step 
algorithm  in  Table  8: 


Table  8:  MATLAB®  Hilbert  Algorithm  [22] 

1.  It  calculates  the  FFT  of  the  input  sequence,  storing  the  result  in  a  vector  x. 

2.  It  creates  a  vector  h  whose  elements  h  ( k)  have  the  values,  where  n  represents 
the  number  of  elements: 

•  1  for  k  =  1,  (n/2 )  +1 

•  2  for  k  =  2,  3,  ...  ,  (n/2 ) 

•  0  for  k  =  (n/2 ) +2,  ...  ,  n 

3.  It  calculates  the  element-wise  product  of  x  and  h. 

4.  It  calculates  the  inverse  FFT  of  the  sequence  obtained  in  step  3  and  returns  the 
first  n  elements  of  the  result. 


Furthermore,  MATLAB®  also  explains  the  uses  of  the  Hilbert  transform  for 
calculations.  The  Hilbert  transform  can  be  used  in  calculating  two  important 
instantaneous  attributes  of  a  signal:  the  amplitude  and  frequency.  The  instantaneous 
amplitude  is  the  amplitude  of  the  complex-valued  Hilbert  transform,  whereas  the 
instantaneous  frequency  is  the  rate  of  change  of  the  instantaneous  phase  angle  with 
respect  to  time  [22] . 

3.3.4.  The  Fourier  Transform 

The  Fourier  transform  is  important  in  many  fields  of  study,  such  as  mathematics, 
engineering,  and  physical  sciences.  Fourier  transforms  are  key  components  in  data 
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processing  and  instruments,  as  well  as  the  cornerstone  of  interferometry  and  aperture 
synthesis  [23].  The  Discrete  Fourier  transform  (DFT),  and  more  specifically  the  FFT,  has 
revolutionized  the  data  processing  and  digital  electronics  industry. 

3.4. 4.1.  Fourier  Transform  Basics 

The  Fourier  transform  is  a  reversible  and  linear  transform  with  a  number  of 
important  properties,  such  as  time  shifting  and  scaling.  For  any  time-domain  function 
fit)  (either  real-  or  complex-valued),  the  Fourier  transform  in  the  frequency  domain  is 
denoted  by  F(n>)  [23]. 

The  forward  transform  of  the  Fourier  transform  is  defined  as: 

/•OO 

Fico)  =  f  f{t)e-JCO,dt  (63) 

J-00 

and  the  inverse  transform  of  the  Fourier  transform  is  represented  as: 

S  r<x> 

fit)  =  —  f  F{co)elwt  dco  .  (64) 

The  complex  exponential  plays  a  key  role  in  the  Fourier  transform. 

A  complex  exponential  is  defined  as  a  complex  number  consisting  of  sinusoids. 
Euler’s  formula  embodies  the  relationship: 

e,e  =  cos(<9)  +  j  sin(<9) .  (65) 

Because  of  the  fact  that  complex  exponentials  are  complete  and  orthogonal,  the 
Fourier  transform  can  represent  any  piecewise  continuous  function  and  to  minimize  the 
error  between  the  function  and  its  transform  representation  [23]. 
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3.4. 4. 2  The  Discrete  Fourier  Transform 


While  the  continuous  Fourier  transform  converts  a  time-domain  signal  of  a 
continuous  spectrum  of  an  infinite  number  of  sinusoids,  the  DFT  is  composed  of  a  finite 
number  of  sinusoids.  The  DFT  is  extremely  important  in  the  area  of  frequency  analysis, 
primarily  due  to  how  it  approaches  a  discrete  signal  in  the  time-domain  and  transforms 
the  signal  into  its  frequency  domain  representation  [24] .  The  DFT  of  a  signal  of  sampled 
data  points,  xn ,  is  defined  by: 


X,  =yN~lx  e~jMIN 

k  ^^Jn=0  n 


and  its  inverse  by: 


(66) 


*n=-YN~lX 

n  AT  /-^k=0 


N  1  y'  jjj'lnnklN 


(67) 


TV  *—‘k=0 

A  DFT  of  N  -point  time  series  results  in  an  N  -point  frequency  spectrum.  When 
the  input  signal  is  real-valued,  the  DFT  contains  an  even  real-part  and  an  odd  imaginary- 
part  of  the  spectrum.  Therefore,  the  “negative”  Fourier  transform  frequencies  provide  no 
new  information,  leading  to  the  conclusion  that  no  information  is  created  nor  destroyed 
by  the  DFT  [23]. 


3.4. 4. 3.  The  Fast-Fourier  Transform 

The  FFT  is  a  faster  version  of  the  DFT.  While  the  FFT  performs  the  same  task  as 
the  DFT,  it  utilizes  algorithms  to  accomplish  the  objectives  in  less  time.  The  discrete 
property  and  speed  of  the  FFT  allows  the  use  of  MATLAB®  to  analyze  signals  [24] . 
MATLAB®’s  f  ft  ( )  function  is  effective  for  computing  the  DFT  of  the  input  signal. 

The  FFT  functions  in  MATLAB®  are  based  on  the  MATLAB®  library,  Fastest 
Fourier  Transform  in  the  West  (FFTW),  is  used  to  increase  the  Fourier  transform  speed. 
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To  compute  the  DFT  of  a  signal,  the  FFTW  library  decomposes  the  problem  using  the 
Cooley-Tukey  algorithm  [22],  resulting  in  the  values  for  the  DFT  of  the  given  signal. 
The  fft  (X)  outputs  the  DFT  of  the  input  data  in  the  form  of  a  vector.  The  output  vector 
is  computed  using  the  FFT  algorithm,  through  selection  using  the  FFTW  library  [22].  To 
further  understand  the  FFT,  the  steps  of  the  MATLAB®  function  fft  ()  can  be  helpful. 
The  function  y=  fft  (X)  implements  the  transform  for  a  vector  of  length  N  as  detailed 
in  below: 


xm=y^mW-'»-n 


(68) 


with 


6)N=e(~2*j),N, 


representing  the  N  th  root  of  unity  [22] . 
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3.5.  Implementation  of  the  Algorithms 

In  the  previous  section,  the  components  and  algorithms  employed  were  explained. 
After  an  understanding  of  the  components  was  gained,  solutions  for  satisfying  the 
objectives  of  this  research  were  implemented  as  code  in  MATLAB®.  Before  the  final 
code  was  implemented  and  used  to  generate  the  results,  there  were  failed  coding  attempts 
explored.  These  attempts  were  detailed  in  the  following  sub- section,  followed  by  the 
final  BEMD  coding  written  for  the  final  results  generated,  and  finishing  the  section  with 
a  detailed  explanation  of  the  implemented  BEMD  code. 
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3.5.1.  Unsuccessful  Coding  Attempts 

There  were  two  unsuccessful  coding  attempts  to  create  the  graphs  to  compare  in 
the  data  analysis  and  results  section.  The  first  attempt  was  to  compare  the  FFT  with  the 
HHT  directly.  Yet,  upon  further  research,  it  was  discovered  that  the  FFT  exists  in  the 
frequency  domain,  whereas  the  HHT  exists  in  the  time  domain.  Therefore,  the  direct 
comparison  attempt  led  to  a  dead  end  due  to  the  inability  to  compare  the  graphs  to  each 
other  since  they  existed  in  different  domains. 

The  next  attempt  was  to  change  the  FFT  into  the  time  domain  using  the  inverse 
FFT  to  compare  directly  with  the  HHT  in  the  same  domain.  The  inverse  FFT  graphs 
presented  an  incorrect  interpretation  of  the  FFT  signal  when  transformed  into  the  time 
domain.  So,  again,  this  second  attempt  was  unsuccessful  in  generating  the  graphs  and 
comparisons  needed  to  successfully  create  the  graphs  for  data  analysis. 

3.5.2.  Bivariate  Empirical  Mode  Decomposition  Explanation 

The  first  component  implemented  in  the  code  was  the  complex  extension  of  the 
EMD  algorithm  BEMD  [11],  described  in  section  3.3.2,  to  calculate  the  complex  IMFs. 
Dr.  Patrick  Flandrin,  one  of  the  primary  authors,  worked  to  create  his  own  code  for 
implementing  the  method  proposed  by  Rilling  et  al.  Through  reading  the  BEMD  paper 
and  proceeding  to  the  website  listed  in  the  paper  (http://perso.ens- 
lyon.fr/patrick.flandrin),  coding  for  the  BEMD  algorithm  was  available  in  the  emd.zip 
file.  Dr.  Flandrin’s  original  BEMD  example  was  detailed  in  Appendix  A. 

3.5.3.  Final  CEMD  Algorithm  Implementation 

Even  though  code  was  written  for  the  BEMD  method,  a  couple  changes  were 
made  to  satisfy  the  objectives  for  this  research.  One  of  the  first  changes  was  to  plot  the 
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magnitude  and  phase  of  the  signal,  rather  than  the  real  and  imaginary  components.  The 
input  complex  RCS  data  set  required  a  change  of  domain  to  produce  useful  phase  and 
magnitude  data  for  graphical  comparison.  Another  modification  was  made  to  present  two 
subplots  on  each  figure,  meaning  that  the  magnitude  and  phase  of  the  decomposed  signal 
were  plotted  as  their  independent  figure.  This  modification  allowed  the  computer  to  run 
without  using  much  memory  with  respect  to  the  computer  used  to  run  the  code,  as  well  as 
allowed  for  a  quicker  side-by-side  comparison  of  the  magnitude  and  phase  plots  of  both 
before  and  after  the  application  of  the  Hilbert  transform.  A  final  change  was  the  labeling 
of  the  title  and  axes,  since  they  need  to  fit  the  real-world  signal  rather  than  the  example 
provided  by  Dr.  Flandrin.  The  modification  of  the  code  written  by  Dr.  Flandrin  is 
presented  in  Appendix  B. 

3. 5. 3.1.  Loading  the  Data  and  Calculation  of  the  IMFs 
Section  B.l,  of  Appendix  B,  represented  the  loading  and  the  remodeling  of 
BEMD  code  to  run  the  data.  Once  the  data  was  loaded,  the  BEMD  algorithm  was  applied 
to  the  remodeled  data.  The  specific  cemdc2  ( )  function  was  presented  in  Appendix  C, 
where  the  definition  of  the  function  was  described  in  detail,  such  as  the  inputs  and  the 
outputs  of  the  function.  After  the  cemdc2  ( )  function  was  employed,  the  function 
outputted  two  separate  values.  The  first,  called  imf ,  was  a  matrix  where  the  IMFs  were 
stored  for  each  radar  return  measurement  of  the  real-world  data  set.  The  second,  called 
nb^iterations,  represented  the  effective  number  of  sifting  iterations  for  each  IMF. 

3.5. 3.2  Plots  of  IMFs:  CEMD  and  HHT 

In  section  B.3  of  Appendix  B,  the  IMFs  that  were  calculated  and  stored  in  section 
B.2  were  plotted.  First,  the  original  signal  was  plotted  with  the  magnitude  on  the  top 


71 


portion  of  the  subplot  and  the  unwrapped  phase  on  the  lower  portion  of  the  subplot.  The 
magnitude  was  plotted  using  the  abs  ( )  function,  whereas  the  phase  was  plotted  using  the 
angle  ( )  function  after  the  value  for  the  angle  was  unwrapped.  The  unwrap  ( )  function 
corrected  the  phase  angles  of  the  input  signal  to  produce  a  smoother  phase  plot  [22]. 
MATLAB®  defines  the  unwrap  ( )  function  as  a  function  that  corrects  the  radian  phase 
angles  in  a  vector  by  adding  multiples  of  +27T  [22].  Once  the  original  signal  was  plotted, 
the  IMFs  and  residue  of  the  decomposed  signal  were  plotted.  The  decomposed  plots 
were  displayed  in  the  same  layout  as  the  original  signal.  The  x-axis  of  the  magnitude  plot 
represented  the  time  scale,  while  the  y-axis  represented  the  RCS  value  in  decibel  per 
square  meter  (dBsm).  This  measurement  represents  the  decibel  measure  for  the  RCS  of  a 
target  relative  to  one  square  meter.  The  x-axis  was  represented  on  the  phase  plot  as  the 
time  scale  as  well,  and  the  y-axis  represented  the  angle  value  of  the  signal  in  radians. 

In  section  B.4  of  Appendix  B,  the  Hilbert  transform  was  applied  to  the  real  and 
imaginary  components  of  the  real-world  data  set.  Rather  than  plot  the  real  and  imaginary 
parts  of  the  data  set,  the  magnitude  and  phase  after  the  application  of  the  Hilbert 
transform  were  plotted.  The  magnitude  plot  represented  the  RCS  value  in  dBsm  on  the  y- 
axis,  with  the  time  scale  represented  on  the  x-axis.  Similarly,  the  phase  plot  represented 
time  on  the  x-axis,  while  the  y-axis  represented  the  angle  of  the  signal  in  radians.  After 
the  Hilbert  transform  of  the  original  signal  was  plotted,  the  Hilbert  transform  was  applied 
to  the  IMFs  and  the  residue  decomposed  from  the  original  data  set,  in  which  the 
magnitude  and  phase  of  the  Hilbert  transform  of  the  IMFs  and  the  residue  were  plotted. 
Because  the  hilbert  ()  function  ignores  the  imaginary  part  of  a  signal,  the  phase  plot 
represents  the  imaginary  component  and  the  magnitude  plot  analyzed  the  real  component 
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of  the  signal.  Similar  to  the  CEMD  plots  in  section  B.3,  the  axes  of  the  Hilbert  transform 
of  the  signal  was  RCS  in  dBsm  with  respect  to  time  for  the  magnitude  plot  and  was  angle 
in  radians  with  respect  to  time  for  the  phase  plot. 

3.5. 3. 3.  Plots  of  IMFs:  FFT  of  CEMD  and  HHT 

The  FFT  of  the  decomposed  signal  before  the  application  of  the  Hilbert  transform 
was  calculated  and  plotted  in  section  B.5  of  Appendix  B.  In  a  similar  manner  to  the 
previous  two  comparisons,  the  original  and  decomposed  signals  were  plotted  for  the  real- 
world  data  set.  Prior  to  generating  the  plots,  the  f ft  ( )  function  was  applied  to  the  real- 
world  data  set.  Then,  the  magnitude  and  phase  were  separately  calculated.  Finally,  the 
magnitude  and  unwrapped  phase  were  plotted.  Since  the  FFT  was  applied  to  the  data  set, 
the  magnitude  graph  represented  frequency  in  GHz  on  the  x-axis  with  the  RCS  value  in 
dBsm  plotted  on  the  y-axis.  Additionally,  the  magnitude  plot  represented  the  absolute 
value  of  the  FFT  of  the  real-world  data  set,  while  using  the  fftshifto  function  applied 
to  the  real-world  data  set  to  shift  the  zero-frequency  component  of  the  data  to  the  center 
of  the  spectrum.  By  using  the  fftshifto  function,  interpretation  of  the  generated  plot 
was  easier.  MATFAB®  states  that  the  fftshift  ()  function  “rearranged  the  outputs  of 
the  [FFT  of  the  data]  by  moving  the  zero-frequency  component  to  the  center  of  the 
array... [and]... is  useful  for  visualizing  a  Fourier  transform  with  the  zero-frequency 
component  in  the  middle  of  the  spectrum”  [22].  The  phase  plot  of  the  FFT  of  the  real- 
world  data  set  was  also  represented  by  frequency  along  the  x-axis  and  the  angle  value  in 
radians  on  the  y-axis. 

After  the  original  signal  was  plotted,  the  plots  of  the  IMFs  and  residue  were 
generated  through  implementation  of  a  for  ()  loop.  The  plots  of  the  IMFs  were  created 
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similar  to  the  original  signal,  with  the  only  exception  of  input  value.  Rather  than  the 
original  data  set,  the  input  signals  were  the  modes  of  the  decomposed  signal  of  the  real- 
world  data  set.  The  magnitude  plot  represented  the  absolute  value  of  the  FFT  of  the 
decomposed  signal,  with  the  fftshift  ()  applied  to  the  resulting  FFT  of  the  IMFs  and 
residue.  For  the  magnitude  plot,  the  x-axis  represented  the  frequency  present  in  the 
modes  of  the  decomposed  signal  and  the  y-axis  represented  the  RCS  value  in  dBsm. 
Additionally,  the  phase  plot  represented  the  angle  of  the  FFT  of  the  modes  of  the 
decomposed  signal.  Frequency  was  represented  along  the  x-axis  and  the  angle  in  radians 
was  represented  along  the  y-axis. 

Section  B.6  detailed  the  process  of  applying  the  FFT  to  the  decomposed  signal 
after  the  application  of  the  Hilbert  transform.  The  first  plots  generated  were  the 
representations  of:  (1)  the  magnitude  of  the  FFT  after  the  Hilbert  transform  application 
and  (2)  the  phase  of  the  FFT  after  the  Hilbert  transform  application.  All  plots  in  this 
section  represented  the  FFT  after  the  application  of  the  Hilbert  transform  of  the  input 
data.  To  generate  the  magnitude  graph,  the  fftshift  ()  of  the  real-valued  data  of  the 
FFT  after  the  application  of  the  Hilbert  transform  was  plotted.  The  x-axis  represented  the 
frequency  component  present  in  the  signal  in  GHz,  while  the  y-axis  represented  the 
amplitude  of  the  RCS  value  in  dBsm.  For  the  phase  plot,  the  unwrapped  angle  of  the 
imaginary-valued  components  FFT  after  the  Hilbert  transform  application  of  the  data  was 
plotted.  The  x-axis  represented  the  amount  of  each  frequency  present  in  the  signal  in 
GHz  and  the  y-axis  represented  the  angle  value  of  the  signal  in  radians. 

After  the  original  signal  of  the  FFT  after  the  application  of  the  Hilbert  transform 
was  generated,  the  magnitude  and  phase  plots  for  the  decomposed  signal  modes  of  the 
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FFT  after  the  Hilbert  transform  application  were  plotted.  The  magnitude  graph  was 
generated  using  the  fftshifto  of  the  real-valued  components  of  the  FFT  after  the 
Hilbert  transform  application  to  the  respective  mode.  On  the  magnitude  graph,  the  x-axis 
represented  the  frequency  present  in  the  signal  and  the  y-axis  represented  the  RCS  value 
in  dBsm.  Additionally,  the  phase  plot  was  generated  through  use  of  the  unwrap  ( )  and 
angle  ( )  functions  from  the  FFT  after  the  Hilbert  transform  application  of  the  imaginary- 
valued  components  of  the  respective  mode.  For  the  phase  plot,  the  x-axis  represented  the 
frequency  present  in  the  signal,  while  the  y-axis  represented  the  angle  of  the  signal  in 
radians. 

3.5. 3.4.  Plots  of  IMFs:  DTI  of  CEMD  and  HHT 

For  the  Doppler-Time-Intensity  (DTI)  plots,  only  the  magnitude  plot  was  of 
interest.  The  DTI  plots  are  represented  on  a  pcolor(),  or  pseudocolor,  plot.  A 
pseudocolor  plot  is  a  rectangular  array  of  cells  with  colors  determined  by  the  input  data 
[22].  The  x-axis  and  y-axis  limits  were  defined  using  the  linspaceO  function  and 
MATLAB®  defines  as  a  function  used  to  “generate  linearly  spaced  vectors”  [22].  To 
graph  a  pcolorO  plot,  there  are  three  dimensions.  The  x-axis  of  the  DTI  was 
represented  by  the  time  value,  the  y-axis  represented  the  frequency  content  in  the  signal 
in  GHz,  and  the  color  axis  represented  the  FFT  intensity  at  each  time  and  frequency  value 
in  the  matrix  spanning  the  x-y  linear  space.  The  third  dimension  was  represented  through 
the  intensity  of  the  color  of  the  absolute  value  of  the  windowed  FFT  both  before  and  after 
the  application  of  the  Hilbert  transform  on  a  logarithmic  scale. 

Section  B.7  in  Appendix  B  presented  the  application  of  the  windowed  FFT  to  the 
original  signal,  IMFs,  and  the  residue  to  generate  the  DTI  plots.  While  a  traditional  FFT 
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would  also  provide  similar  analysis,  the  FFT  employed  in  this  section  of  the  code  was  a 
windowed  FFT  implemented  by  the  user,  rather  than  a  function  embedded  in  the 
MATLAB®  program.  The  windowed  FFT  replaces  the  Fourier  transform's  sinusoidal 
wave  by  the  product  of  a  sinusoid  and  a  window  which  is  localized  in  time.  In  addition, 
windowed  Fourier  Transforms  are  important  in  providing  simultaneous  insight  in  time 
and  frequency  behavior  of  the  functions.  The  first  step  was  to  define  the  step  resolution 
and  the  amount  of  seconds  for  the  resolution  of  the  windowed  FFT.  Next,  the  windowed 
FFT  was  applied  and  the  DTI  plot  of  the  FFT  of  the  original  signal  was  generated. 

Subsequently,  the  windowed  FFT  of  the  decomposed  signal  modes  were 
calculated  and  stored.  Then,  the  FFT  of  the  decomposed  signal  modes  were  graphed  on  a 
pcolor  ( )  plot.  These  two  objectives  were  accomplished  using  an  embedded  for  ( )  loop. 
First,  the  windowed  FFT  value  of  the  IMF  was  calculated  and  the  magnitude  of  the  FFT 
was  graphed.  The  fftshift  ()  of  the  absolute  value  of  the  windowed  FFT  was  graphed 
on  the  pcolor  ( )  plot  using  a  logarithmic  scale.  The  DTI  plot  was  represented  with  time 
on  the  x-axis,  frequency  on  the  y-axis,  and  the  color  axis  representing  the  intensity  of  the 
FFT  of  the  decomposed  signal. 

Section  B.8  of  Appendix  B  represented  the  calculation  of  the  windowed  FFT  of 
the  original  signal  and  the  modes  of  the  decomposed  signal  after  the  application  of  the 
Hilbert  transform,  as  well  as  generated  the  DTI.  The  first  step  was  to  define  the  step 
resolution  and  the  amount  of  seconds  for  the  resolution  of  the  windowed  FFT.  Next,  the 
windowed  FFT  after  the  Hilbert  transform  application  of  the  original  signal  was 
calculated  and  stored.  Once  the  windowed  FFT  was  calculated  and  store,  the  windowed 
FFT  was  plotted. 
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The  windowed  FFT  of  each  mode  of  the  decomposed  signal  was  calculated  and 
stored.  Then,  the  FFT  after  the  Hilbert  transform  application  of  the  decomposed  signal 
graphed  on  a  pcolor()  plot.  These  two  objectives  were  accomplished  using  an 
embedded  for()  loop.  First,  the  Hilbert  transform  of  each  component  of  the 
decomposed  signal  was  calculated  and  the  windowed  FFT  value  after  the  Hilbert 
transform  application  was  calculated.  Next,  the  magnitude  of  the  FFT  after  the  Hilbert 
transform  application  was  plotted.  The  fftshifto  of  the  absolute  value  of  the 
windowed  FFT  after  the  Hilbert  transform  application  was  graphed  on  the  pcolor  ( )  plot 
using  a  logarithmic  scale.  The  DTI  plot  was  represented  with  time  on  the  x-axis, 
frequency  on  the  y-axis,  and  the  color  axis  representing  the  intensity  of  the  FFT  after  the 
Hilbert  transform  application  of  the  decomposed  signal. 

3.6.  Chapter  Summary 

An  exploration  of  the  various  components  necessary  to  implement  the  MATLAB® 
code  and  satisfy  the  objectives  presented  in  Chapter  I  was  performed.  The  applied 
BEMD  algorithm  was  explained  and  relevant  data  collection  was  detailed.  Four  specific 
components  were  investigated  for  employing  the  algorithm  presented  by  Rilling  et  al. 
[11]:  (1)  the  EMD  algorithm,  (2)  the  CEMD  algorithm,  (3)  the  Hilbert  transform  (leading 
to  the  HHT  representation),  and  (4)  the  FFT.  Finally,  algorithm  was  implemented  and 
discussed  with  detailed  MATFAB®  code  as  presented  in  Appendices  A,  B,  and  C. 
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IV.  Results  and  Data  Analysis 


4.1.  Overview 

The  results  and  analysis  for  the  comparisons  of  the  data  generated  in  Chapter  III, 
after  applying  the  CEMD  algorithm  and  the  two  mathematical  transforms,  the  Hilbert 
transform  and  FFT,  are  presented  and  analyzed  in  detail.  The  results  are  organized  into 
four  sections.  The  first  section  compares  the  decomposed  real-world  data  set  with  the 
RCS  signal  of  various  canonical  shapes  collected  in  the  RCS  range  at  AFIT,  often 
encountered  in  real-world  signals.  Next,  the  decomposed  real-world  signal  is  compared 
with  its  Hilbert  transform  of  the  decomposed  real-world  signal.  In  addition,  the 
application  of  the  FFT  to  both  the  decomposed  signal  and  the  application  of  the  Hilbert 
transform  to  the  decomposed  signal  are  compared.  Finally,  the  windowed  FFT  of  the 
decomposed  signal  and  the  Hilbert  transform  application  to  the  decomposed  signal  are 
compared.  In  other  words,  the  Doppler-Time-Intensity  (DTI)  representation  of  the 
signals  were  analyzed  and  compared,  noting  the  differences  and  similarities  in  the 
generated  graphical  representation. 

4.2.  Evaluation  of  CEMD  Method 

To  assess  how  well  the  CEMD  method  performed,  the  decomposed  signal  of  the 
real-world  data  set  was  compared  with  four  common  shapes  inherent  to  dynamic  objects 
similar  to  the  one  represented  by  the  real-world  data  set.  First,  the  graph  of  the  original 
signal  of  a  closed-cap  cylinder  was  compared  to  the  various  IMFs  of  the  decomposed 
real-world  signal.  After  comparing  the  RCS  of  the  closed-cap  cylinder  with  the  IMFs, 
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the  most  similar  IMF  to  the  original  signal  of  the  closed-cap  cylinder  was  IMF  #2,  as 


displayed  in  Figures  8  and  9. 


CEMD  Magnitude  of  Closed-Cap  Cylinder-Original  Signal 
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Figure  8:  CEMD  Magnitude  of  Closed-Cap  Cylinder 


CEMD  Magnitude  of  Real-World  Signal— IMF  #2 
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Figure  9:  CEMD  Magnitude  of  Real-World  Signal— IMF  #2 


The  peaks  in  the  RCS  values  of  the  two  plots  occurred  at  analogous 
fixed  points  in  the  analysis  of  the  complex  time  series.  However,  the  amplitude  of  the 
peak  return  values  differ  by  a  factor  of  20  dB.  Even  though  this  difference  exists,  the 
RCS  value  of  the  original  signal  of  the  closed-cap  cylinder  behaved  most  like  IMF  #2  of 
the  decomposed  real-world  signal  through  comparing  the  plots. 

The  second  comparison  occurred  between  the  RCS  of  a  cavity  cylinder  and  the 
decomposed  real-world  signal.  The  cavity  cylinder  results  were  similar  to  the  closed-cap 
cylinder  RCS  return  plot.  Additionally,  the  peaks  of  the  cavity  cylinder  and  IMF  #2 
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occur  at  approximately  the  same  instances  in  time  on  the  graph,  shown  in  Figures  10  and 


11. 


CEMD  Magnitude  of  Cavity  Cylinder-Original  Signal 
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Figure  10:  CEMD  Magnitude  of  Closed-Cap  Cylinder 


CEMD  Magnitude  of  Real-World  Signal— IMF  #2 


Figure  11:  CEMD  Magnitude  of  Real-World  Signal— IMF  #2 


The  amplitude  of  the  peaks  was  also  different  by  a  factor  of  20  dB;  although  there 
was  a  difference,  the  cavity  cylinder  behaved  most  like  IMF  #2  of  the  decomposed  real- 
world  signal. 

Third,  the  RCS  plot  of  the  cone-sphere  was  compared  with  the  IMFs  of  the  real 
world  signal.  Yet,  unlike  the  first  two  shape  comparisons,  there  was  no  IMF  of  the 
decomposed  real-world  signal  that  behaved  most  like  the  RCS  plot  of  the  cone-sphere 
signal. 
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Finally,  the  RCS  plot  of  the  original  signal  of  the  dihedral  corner  reflector  was 
compared  to  the  IMFs  of  the  decomposed  real-world  signal.  Comparing  these  two  plots 
side-by-side,  the  dihedral  corner  reflector  behaved  most  like  IMF  #3,  with  a  magnitude 
difference  by  a  factor  of  seven.  One  possible  explanation  for  the  difference  in  the 
magnitude  factors  is  the  amount  of  space  between  the  dynamic  object  and  the  radar 
measuring  the  RCS  value.  Finally,  the  peaks  of  each  plot  occurred  at  approximately  the 
same  proportional  time  interval,  and  the  smaller  radar  returns  also  occurred  at  the 
approximate  same  time  intervals. 


81 


CEMD  Magnitude  of  Real-World  Signal— IMF  #3 
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Figure  14:  CEMD  Magnitude  of  Real-World  Signal — IMF  #3 


Overall,  the  comparisons  discussed  reinforced  the  fact  that  the  CEMD  method 
effectively  decomposed  the  real-world  data  set  provided  for  this  evaluation  into  its  IMFs. 
Therefore,  the  CEMD  method  could  be  a  viable  tool  for  analysis  of  RCS  signals  in  the 
future. 


4.3.  Analysis  of  Decomposed  Data  Before  and  After  the  Hilbert  Transform 

The  second  analysis  was  the  comparison  of  the  graphical  results  generated  of  the 
original  signal  and  the  Hilbert  transform  application  to  the  original  signal  using  the 
CEMD  algorithm.  There  was  not  any  quantitative  data  to  analyze;  rather,  the  graphs 
were  qualitatively  compared  and  discussed. 
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4.3.1.  Original  Signal  Plot  Analysis 


Firstly,  the  magnitude  of  the  original  signal  and  the  Hilbert  transform  application 


to  the  original  signal  were  displayed  in  Figures  15  and  16  and  represented  by  the  graph  of 


RCS  value  in  dBsm  with  respect  to  time. 


CEMD  Magnitude  of  Real-World  Signal-Original  Signal 
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Figure  15:  CEMD  Magnitude  of  Original  Signal 


50 

40 

30 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal-Original  Signal 
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Figure  16:  HHT  Magnitude  of  Original  Signal 


The  graphs  both  displayed  five  large  returns,  where  the  original  signal  had  a 
slightly  higher  return  at  each  spike  than  after  the  application  of  the  Hilbert  transform. 
Other  than  the  differences  in  the  peaks,  there  was  no  noticeable  difference  in  the  original 
signal  when  compared  with  the  signal  after  the  application  of  the  Hilbert  transform. 

However,  the  phase  plots,  shown  in  Figures  17  and  18,  displayed  a  much  larger 
difference.  Namely,  the  scale  on  the  y-axis  displayed  a  large  difference  in  amplitude,  as 
well  as  the  smoothness  of  the  curve.  The  plots  needed  to  be  viewed  on  the  same  scale  to 
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determine  if  there  was  a  difference  in  the  curve  because  the  plots  were  on  such  different 


scales. 


CEMD  Phase  of  Real-World  Signal-Original  Signal 


Figure  17:  CEMD  Phase  of  Original  Signal 

Hilbert-Huang  Transform  Phase  of  Real-World  Signal-Original  Signal 


Figure  18:  HHT  Phase  of  Original  Signal 


However,  these  differences  and  possible  causes  fell  outside  the  scope  of  this 
evaluation  and  would  require  more  research  to  determine  whether  the  phase  portion  of 
the  CEMD  algorithm  provided  significant  analysis. 

Comparatively  to  the  magnitude  plots  of  the  original  signal,  there  were  not 
noticeable  differences  between  the  before  and  after  the  application  of  the  Hilbert 
transform.  Consequently,  more  information  and  noticeable  changes  might  occur  after 
analysis  of  specific  IMFs. 
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4.3.2.  IMF  Plot  Analysis 

The  plots  of  IMF  #1  showed  only  a  small  difference  between  the  two  plots  both 
before  and  after  the  application  of  the  Hilbert  transform,  in  Figures  19  and  20.  Therefore, 
no  helpful  information  could  be  determined  by  comparison  of  these  two  plots. 


CEMD  Magnitude  of  Real-World  Signal— IMF  #1 
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Figure  19:  EMD  Magnitude  of  IMF  #1 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #1 


20 

»  15 

CD 

w  10 

o 

*  5 

_ L 

1 . 1 . - . - 

2000  4000  6000  8000  10000  12000  14000 

Time 

Figure  20:  HHT  Magnitude  of  IMF  #1 


By  comparing  the  two  IMF  #2  plots  in  Figures  21  and  22,  there  was  a  slight 
difference  in  the  intensity  of  each  return  value  with  respect  to  time.  On  the  graph  of  the 
IMF  #2  signal  before  the  application  of  the  Hilbert  transform,  there  was  approximately 
double  the  intensity  in  time  when  compared  with  the  signal  after  the  application  of  the 
Hilbert  transform.  Although  there  was  a  slight  magnitude  difference,  the  primary  and 
most  noticeable  difference  was  the  intensity  at  each  time  value  before  the  application  of 
the  Hilbert  transform  signal.  After  application  of  the  Hilbert  transform,  the  negative 
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frequency  values  were  set  to  equal  zero  and  could  account  for  the  less  intense  plot  value 


after  the  Hilbert  transform  application.  For  a  more  detailed  representation  of  what  the 


Hilbert  transform  does  to  an  input  signal,  please  refer  to  Chapter  III,  section  3.3.3. 
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Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #2 


Beginning  with  IMF  #5,  displayed  in  Figures  23  and  24,  the  curve  after 
application  of  the  Hilbert  transform  presented  a  more  smooth  representation  than 
previous  IMFs,  as  well  as  displayed  a  more  precise  value  representation  of  the  original 
signal.  IMF  #5  provided  a  cleaner  representation  of  the  decomposed  real-world  signal, 
even  though  the  amplitude  scale  was  slightly  different  between  two  plots. 
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Figure  21:  CEMD  Magnitude  of  IMF  #2 
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Figure  23:  CEMD  Magnitude  of  IMF  #5 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #5 


From  comparing  the  plots  of  IMF  #6,  in  Figures  25  and  26,  this  IMF  allowed  for 
direct  comparison  and  showed  that,  after  application  of  the  Hilbert  transform,  the  plot 
was  more  fluid  than  before  the  application  of  the  Hilbert  transform.  The  same 
characteristics  were  present  in  both  plots  and,  through  application  of  the  Hilbert 
transform,  the  plot  of  IMF  #6  of  the  real-world  signal  was  more  focused  and  displayed  a 
cleaner  representation.  The  “cleaned-up”  version  makes  analysis  of  the  time  value 
representation  in  RCS  value  easier  and  may  have  provided  more  accurate  analysis  of  the 
decomposed  signal  after  the  application  of  the  Hilbert  transform. 
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CEMD  Magnitude  of  Real-World  Signal— IMF  #6 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #6 


IMFs  #7  through  #11  provided  additional  support  for  the  conclusions  and 
characteristics  discussed  about  IMF  #6. 

The  decomposition  displayed  in  Figures  27  through  36  showed  the  sinusoidal 
nature  of  the  CEMD  method;  however,  the  sinusoid  became  more  noticeable  after 
application  of  the  Hilbert  transform.  Before  the  application  of  the  Hilbert  transform,  the 
decomposition  of  the  signal  appeared  noisier  and  possibly  contained  extraneous  signal 
values  of  the  decomposed  signal  representation.  Based  on  the  aforementioned 
conclusions,  it  is  hypothesized  in  later  plot  representations  that  the  removal  of  negative 
frequencies  might  prove  be  either  helpful  or  harmful,  especially  when  used  for  analysis  in 
the  intelligence  community. 
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CEMD  Magnitude  of  Real-World  Signal— IMF  #7 


Figure  27:  CEMD  Magnitude  of  IMF  #7 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #7 


CEMD  Magnitude  of  Real-World  Signal— IMF  #8 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #8 


Figure  30:  HHT  Magnitude  of  IMF  #8 
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CEMD  Magnitude  of  Real-World  Signal— IMF  #9 


Figure  31:  CEMD  Magnitude  of  IMF  #9 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #9 


Figure  32:  HHT  Magnitude  of  IMF  #9 


CEMD  Magnitude  of  Real-World  Signal— IMF  #10 


Figure  33:  CEMD  Magnitude  of  IMF  #10 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #10 


Figure  34:  HHT  Magnitude  of  IMF  #10 
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CEMD  Magnitude  of  Real-World  Signal-IMF  #11 


Figure  35:  CEMD  Magnitude  of  IMF  #11 


Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal— IMF  #11 


4.4.  Analysis  of  the  FFT  Plots  Before  and  After  the  Hilbert  Transform  Application 

Traditionally,  the  FFT  has  been  the  main  signal  analysis  tool  for  analyzing  RCS 
data.  The  application  of  the  FFT  allowed  the  researcher  to  determine  at  what  frequency 
values  the  signal  was  most  dominant,  as  well  as  represented  the  signal  for  analysis  in  the 
frequency  domain.  For  each  of  the  graphs  generated  to  represent  the  magnitude  of  the 
FFT  in  the  frequency  domain,  the  frequencies  represented  by  7685  Hz  and  below  were 
the  negative  frequencies,  while  the  positive  frequencies  were  between  7686  Hz  and 
15370  Hz  on  the  frequency  axis.  In  addition  to  the  magnitude  plot  being  generated,  the 
phase  of  the  signal  was  represented  by  angle  in  radians  with  respect  to  the  frequency  of 
the  signal. 
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4.4.1.  Original  Signal  FFT  Plot  Comparisons 

The  first  plots  compared  the  effect  of  the  Hilbert  transform  applied  to  the  real- 
world  signal.  The  magnitude  plots  of  the  original  signal  before  and  after  applying  the 
Hilbert  transform,  prior  to  decomposition  using  the  CEMD  algorithm,  are  represented 
below  in  Figures  37  and  38. 


FFT  Magnitude  of  Real-World  Signal-Original  Signal 
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Figure  37:  FFT  of  CEMD  Magnitude  of  Original  Signal 


HHT  of  FFT  Magnitude  of  Real-World  Signal-Original  Signal 
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Figure  38:  FFT  of  HHT  Magnitude  of  Original  Signal 


First,  the  comparison  of  the  magnitude  plot  of  the  FFT  before  and  after  the  Hilbert 
transform  application  displayed  a  characteristic  inherent  to  the  Hilbert  transform.  The 
characteristic  inherent  to  the  Hilbert  transform  resulted  in  the  negative  frequencies  of  the 
signal  being  removed,  leading  the  negative  frequency  values  to  equal  zero.  The  FFT  of 
the  signal  before  the  Hilbert  transform  application  showed  the  plot  contained  both 
negative  and  positive  frequencies,  with  the  signal  appearing  approximately  symmetric 
along  the  frequency  axis  with  respect  to  its  RCS  value  in  dBsm.  The  removal  of  the 
negative  frequencies  on  the  plot  was  not  the  only  difference  present  in  the  generated 
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graph.  The  RCS  return  values  that  were  present  in  the  positive  frequencies  after  the 
application  of  the  Hilbert  transform  and  the  graph  had  a  larger  magnitude  value  by  a 
factor  of  approximately  two  as  compared  with  the  RCS  return  values  before  the 
application  of  the  Hilbert  transform.  The  difference  in  amplitude  of  the  RCS  return  value 
led  to  the  possibility  that  the  negative  frequency  returns  were  preserved  in  the  graph  and 
resulted  in  the  larger  magnitude  present  on  the  positive  frequency  portion  of  the  plot. 
One  possible  explanation  is  based  on  the  conservation  of  energy  principle,  which  states 
that  the  amount  of  energy  put  into  the  signal  must  equal  that  exiting  after  application  of 
the  Hilbert  transform. 

Additionally,  the  phase  plots  of  the  original  signal  of  the  FFT  plot  both  before  and 
after  the  Hilbert  transform  application  did  not  appear  to  provide  any  new  or  possibly 
significant  information  due  to  them  displaying  approximately  equal  angle  measurements, 
presented  in  Figures  39  and  40. 


FFT  Phase  of  Real-World  Signal-Original  Signal 


Figure  39:  FFT  of  CEMD  Phase  of  Original  Signal 
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HHT  of  FFT  Phase  of  Real-World  Signal-Original  Signal 


Figure  40:  FFT  of  FIHT  Phase  of  Original  Signal 


4.4.2.  IMF  FFT  Plot  Comparisons 

In  a  similar  fashion  to  the  plots  of  the  original  signal,  comparisons  between  the 
IMFs  of  the  FFT  after  application  of  the  Hilbert  transform  of  the  decomposed  signal  yield 
different  results  and  noted  graphical  differences.  The  plots  representing  the  IMFs  of  the 
FFT  were  detailed  below  in  Figures  41  through  50. 

The  main  difference  between  the  IMF  plots  of  the  FFT  plot  before  and  after  the 
Hilbert  transform  was  the  disappearance  of  the  negative  frequencies  after  the  application 
of  the  Hilbert  transform,  whereas  before  the  application  of  the  Hilbert  transform,  the  full 
frequency  spectrum  displayed  the  radar  return  signal.  Removal  of  the  negative 
frequencies  of  the  signal  after  application  of  the  Hilbert  transform  is  a  property  of  the 
FFT  of  the  Hilbert  transform  and  discussed  in  Section  3.3.3.  In  addition  to  the  difference 
along  the  frequency  axis,  the  magnitude  of  the  radar  return  varied  between  the  graphs 
generated  before  and  after  the  application  of  the  Hilbert  transform. 
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FFT  Magnitude  of  Real-World  Signal— IMF  #1 
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Figure  41:  FFT  of  CEMD  Magnitude  of  IMF  #1 


HHT  of  FFT  Magnitude  of  Real-World  Signal-IMF  #  1 


FFT  Magnitude  of  Real-World  Signal-IMF  #2 
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Figure  44:  FFT  of  HHT  Magnitude  of  IMF  #2 
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FFT  Magnitude  of  Real-World  Signal— IMF  #3 


HHT  of  FFT  Magnitude  of  Real-World  Signal-IMF  #3 
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Figure  46:  FFT  of  HHT  Magnitude  of  IMF  #3 
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FFT  Magnitude  of  Real-World  Signal-IMF  #4 
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Figure  47:  FFT  of  CEMD  Magnitude  of  IMF  #4 

HHT  of  FFT  Magnitude  of  Real-World  Signal-IMF  #4 
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FFT  Magnitude  of  Real-World  Signal— IMF  #5 
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Figure  49:  FFT  of  CEMD  Magnitude  of  IMF  #5 

HHT  of  FFT  Magnitude  of  Real-World  Signal-IMF  #5 


Prior  to  the  Hilbert  transform  application,  IMF  plots  were  consistent  with  a  high- 
end  limit  of  the  complex  RCS  return  value  between  40  dB  and  approximately  250  dB 
from  IMF  #1  to  IMF  #11.  Also,  the  magnitude  plots  after  the  Hilbert  transform  plots 
ranged  with  a  high-end  limit  from  75  dB  to  310  dB,  with  IMFs  #6  through  #11  are 
detailed  in  Figures  51  through  62.  The  difference  in  magnitude  was  indicative  of  the 
conclusion  that  the  HHT  was  not  a  useful  analysis  tool  in  comparison  with  the  FFT  for 
the  real-world  data  set  provided  for  this  project.  A  deeper  understanding  of  exactly  what 
the  Hilbert  transform  application  yields  from  a  similar  complex  RCS  data  set  is  necessary 
for  determining  the  usability  of  the  HHT  as  a  tool  for  complex  RCS  data  used  in  the 
intelligence  community,  when  analyzed  in  the  frequency  domain. 
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Figure  51:  FFT  of  CEMD  Magnitude  of  IMF  #6 
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Figure  52:  FFT  of  HHT  Magnitude  of  IMF  #6 
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Figure  54:  FFT  of  HHT  Magnitude  of  IMF  #7 
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FFT  Magnitude  of  Real-World  Signal— IMF  #8 
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Figure  55:  FFT  of  CEMD  Magnitude  of  IMF  #8 
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Figure  56:  FFT  of  HHT  Magnitude  of  IMF  #8 
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Figure  57:  FFT  of  CEMD  Magnitude  of  IMF  #9 
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Figure  58:  FFT  of  HHT  Magnitude  of  IMF  #9 
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FFT  Magnitude  of  Real-World  Signal— IMF  #10 
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Figure  59:  FFT  of  CEMD  Magnitude  of  IMF  #10 
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Figure  60:  FFT  of  HHT  Magnitude  of  IMF  #10 
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Figure  61:  FFT  of  CEMD  Magnitude  of  IMF  #11 
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Figure  62:  FFT  of  HHT  Magnitude  of  IMF  #11 


The  phase  plots  of  the  IMFs  of  the  FFT  plot  of  before  and  after  the  Hilbert 
transform  application  provided  limited  new  or  possibly  significant  data,  although  the 
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plots  displayed  a  significant  difference  between  the  two  mathematical  transforms,  with 


the  IMFs  represented  in  Figures  63  through  70. 
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Figure  63:  FFT  of  CEMD  Phase  of  IMF  #1 
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Figure  64:  FFT  of  HHT  Phase  of  IMF  #1 
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Figure  65:  FFT  of  CEMD  Phase  of  IMF  #5 
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Figure  66:  FFT  of  HHT  Phase  of  IMF  #5 
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FFT  Phase  of  Real-World  Signal-IMF  #7 


Figure  67:  FFT  of  CEMD  Phase  of  IMF  #7 
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Figure  68:  FFT  of  HHT  Phase  of  IMF  #7 
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Figure  69:  FFT  of  CEMD  Phase  of  IMF  #11 
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Figure  70:  FFT  of  HHT  Phase  of  IMF  #11 


A  distinct  difference  between  these  phase  plots  both  before  and  after  the 
application  of  the  Hilbert  transform  existed;  however,  the  understanding  needed  to 
determine  the  significance  or  lack  thereof  of  these  plots  is  outside  the  scope  of  this 
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inquiry.  Further  understanding  of  the  phase  plots  and  the  interpretation  of  these  plots  is 
necessary  to  determine  new  or  significant  information. 

4.5.  Analysis  of  Doppler-Time-Intensity  Plots  Prior  to  and  After  the  Hilbert 
Transform 

4.5.1.  Doppler-Time-Intensity  (DTI)  Plots 

For  intelligence  analysts,  the  DTI  plot  is  of  interest  for  understanding  what 
occurred  at  each  frequency  and  time  interval,  as  well  as  displaying  the  amount  of  energy 
being  transmitted  at  that  moment.  This  plot  is  created  using  a  windowed  Fourier 
transform  and  is  an  accepted  tool  in  the  analyst  community.  More  information 
concerning  DTI  plots  are  in  Section  3. 5. 3. 4. 

4.5.2.  Original  Signal  DTI  Plot  Comparison 

Immediately,  a  major  difference  was  noticed  in  the  plots  between  the  original  DTI 
plots  of  before  and  after  the  Hilbert  transform  application,  in  Figures  71  and  72.  Prior  to 
the  application  of  the  Hilbert  transform,  the  signal  displayed  a  symmetric  signal  along  the 
frequency  axis,  whereas  the  negative  frequency  values  after  the  Hilbert  transform 
application  were  no  longer  symmetric  and  did  not  allow  for  analysis  capabilities  in  the 
known  intelligence  community.  Representation  of  the  HHT  as  a  DTI  plot  yielded  reason 
to  believe  the  application  of  the  Hilbert  transform  did  not  allow  for  the  proper  analysis 
needed  for  RCS  data. 
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Figure  71:  DTI  Magnitude  of  Original  Signal 
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Figure  72:  DTI  HHT  Magnitude  of  Original  Signal 


The  time  axis  from  0  to  50  displayed  a  more  intense  return  on  the  DTI  plot  after 
the  Hilbert  transform  application,  according  to  the  color  representation.  When  compared 
with  the  DTI  plot  of  before  the  application  of  the  Hilbert  transform  signal,  the  plot 
representation  after  the  application  of  the  Hilbert  transform  provided  a  possibly  incorrect 
interpretation  of  the  signal.  Thus,  the  only  useable  portion  of  the  DTI  plot  after  the 
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Hilbert  transform  application  of  the  pre-decomposed  signal  was  the  positive  frequency 
between  time  50  and  time  100.  Nevertheless,  before  the  decomposition  of  the  signal 
using  the  CEMD  algorithm,  the  application  of  the  Hilbert  transform  did  not  provide  a 
DTI  plot  that  yielded  useful  results. 

4.5.3.  IMF  DTI  Plot  Comparisons 

When  comparing  the  two  plots  of  IMF  #1,  represented  in  Figures  73  and  74,  the 
same  characteristics  and  differences  present  in  the  comparison  of  the  original  signal  were 
also  present  in  this  specific  IMF.  Again,  after  the  application  of  the  Hilbert  transform, 
the  negative  frequencies  were  non-existent  and  little  useful  information  was  displayed  on 
the  plot  of  IMF  #1.  The  large  returns,  normally  indicative  of  broadside  flashes  in  RCS 
data,  were  shifted  in  the  DTI  plot  after  the  Hilbert  transform  application  and  potentially 
resulted  from  the  fftshift  ()  or  the  fft  ()  functions  in  MATFAB®.  This  shifting  was 
not  noticeable  until  IMF  #1  and  many  techniques  were  employed  to  attempt  to  fix  this 
shifting;  however,  the  fftshift  ()  and  the  ifftshifto  functions  were  not  helpful  with 
shifting  the  time  scale,  where  the  problem  appeared  to  occur.  The  time  shifting  occurred 
after  the  application  of  the  Hilbert  transform  presented  another  reason  why  the  Hilbert 
transform  application  to  the  real-world  signal  did  not  produce  a  more  enhanced  fidelity  or 
a  better  analysis  tool  for  RCS  data. 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #1 
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Figure  73:  DTI  Magnitude  of  IMF  #1 
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Figure  74:  DTI  HHT  Magnitude  of  IMF  #1 


From  the  DTI  plots  for  IMF  #2,  in  Figures  75  and  76,  both  before  and  after  the 
Hilbert  transform  application,  the  plots  displayed  differences  detectable  between  the 
before  and  after  the  application  of  the  Hilbert  transform  plots.  The  negative  frequencies 
were  unusable,  even  though  the  plot  started  looking  more  like  a  mirror  image  of  the 
positive  frequencies.  The  characteristic  of  growing  intensity  on  the  negative  frequency 
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portion  of  the  graph  may  represent  the  higher,  more  intense  mode  being  sifted  out  of  the 
decomposed  signal  on  the  plot.  The  time  axis  was  still  shifted  after  the  application  of  the 
Hilbert  transform.  The  shifting  might  be  caused  by  a  function  inherent  to  the  Hilbert 
transform,  which  displayed  another  reason  that  the  application  of  the  Hilbert  transform 
may  not  provide  new  or  enhanced  information  than  the  application  of  the  FFT. 
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Figure  75:  DTI  Magnitude  of  IMF  #2 


DTI  Plot-HHT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #2 
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Figure  76:  DTI  HHT  Magnitude  of  IMF  #2 
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By  comparing  IMFs  #5  through  #8,  the  negative  frequencies  after  the  Hilbert 
transform  application  on  the  plots  became  less  intense,  represented  by  Figures  77  through 
84;  however,  the  higher  and  more  intense  RCS  values  were  being  sifted  out  using  the 
CEMD  algorithm,  and  provided  positive  frequencies  mirrored  as  negative  frequencies. 
When  comparing  the  DTI  plots  before  the  Hilbert  transform  and  the  signal  after  the 
application  of  the  Hilbert  transform  characteristics  present  prior  to  the  application  of  the 
Hilbert  transform  were  not  represented.  Therefore,  if  the  Hilbert  transform  is  applied  to 
the  data,  data  is  lost,  leading  to  the  conclusion  that  the  analysis  would  not  represent  the 
true  signal.  By  analyzing  these  six  IMFs,  the  Hilbert  transform  application  appeared  to 
ineffectively  represent  the  data  decomposed  using  the  CEMD  algorithm. 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #5 
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Figure  77:  DTI  Magnitude  of  IMF  #5 
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Figure  78:  DTI  HHT  Magnitude  of  IMF  #5 
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Figure  79:  DTI  Magnitude  of  IMF  #6 


DTI  Plot-HHT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #6 


Time  (s) 

Figure  80:  DTI  HHT  Magnitude  of  IMF  #6 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #7 
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Figure  81:  DTI  Magnitude  of  IMF  #7 
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Figure  82:  DTI  HHT  Magnitude  of  IMF  #7 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #8 
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Figure  83:  DTI  Magnitude  of  IMF  #8 


Figure  84:  DTI  HHT  Magnitude  of  IMF  #8 


IMFs  #9  through  #1 1  display  major  differences  in  the  DTI  plots  before  and  after 
the  Hilbert  transform  application,  shown  in  Figures  85  through  90. 
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First,  the  shift  of  the  time  axis  was  noticeable  in  IMF  #9  due  to  the  vertical  line 


present  at  time  equal  to  50.  Secondly,  the  signal  displayed  before  the  Hilbert  transform 
was  more  focused  and  displayed  all  data  points,  whereas  the  graph  after  the  application  of 
the  Hilbert  transform  was  smoother  and  overlooked  the  individual  spikes  from  the  signal 
decomposition.  Finally,  it  was  concluded  that  the  signal  after  the  Hilbert  transform 
application  was  not  representative  of  the  originally  decomposed  signal. 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #9 
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Figure  85:  DTI  Magnitude  of  IMF  #9 


DTI  Plot- FI  FIT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #9 
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Figure  86:  DTI  HHT  Magnitude  of  IMF  #9 
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Figure  87:  DTI  Magnitude  of  IMF  #10 
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Figure  88:  DTI  HHT  Magnitude  of  IMF  #10 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #11 
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Figure  89:  DTI  Magnitude  of  IMF  #11 

DTI  Plot-HHT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #  1 1 
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Figure  90:  DTI  HHT  Magnitude  of  IMF  #11 


4.5.4.  DTI  Plot  Comparison  Conclusions 

By  comparing  the  DTI  plots  both  before  and  after  the  application  of  the  Hilbert 
transform  to  the  original  signal  and  decomposed  data,  it  was  concluded  the  application  of 
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the  HHT  did  not  provide  new  or  significant  enhancements  to  the  analysis  of  the  data. 
After  the  detailed  comparison  of  each  IMF,  the  Hilbert  transform  application  seemed  to 
delete  approximately  half  of  the  original  frequency  data,  as  well  as  incorrectly 
represented  the  signal  along  the  time  scale  axis.  It  is  hypothesized  that  the  HHT  might 
provide  a  more  enhanced  fidelity  and  new  or  useful  information  of  the  input  RCS  signal. 
However,  this  application  of  the  Hilbert  transform  did  not  provide  a  more  enhanced 
fidelity  than  the  application  of  the  FFT.  Additionally,  the  Hilbert  transform  removed 
important  data  and  did  not  display  the  correct  amount  of  fidelity  compared  to  the 
application  of  the  FFT.  It  cannot  be  clearly  determined  whether  the  HHT  was  a  useful 
signal  analysis  tool  for  RCS  data.  Furthermore,  it  could  not  be  determined  that,  when 
compared  with  the  application  of  the  FFT,  the  Hilbert  transform  was  a  better 
representation  of  the  data  and  was  more  a  hindrance  than  helpful  tool  for  the  analysis  of 
the  DTI  plots. 

4.6.  Chapter  Summary 

First,  the  magnitude  of  the  decomposed  signal  was  compared  with  various 
canonical  RCS  graphs  to  assess  the  CEMD  algorithm  in  decomposing  the  real-world  data 
set.  Next,  the  magnitude  of  the  decomposed  signal  was  compared  both  before  and  after 
the  application  of  the  Hilbert  transform.  Similarly,  the  magnitude  of  the  FFT  of  the 
decomposed  signal  was  compared  both  before  and  after  the  application  of  the  Hilbert 
transform.  Finally,  the  generated  magnitude  DTI  plots  of  the  decomposed  signal  before 
and  after  the  application  of  the  Hilbert  transform  were  compared.  In  addition,  the  phase 
plots  for  the  decomposed  signal  and  the  FFT  of  the  decompose  signal  were  compared 
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with  their  respective  graphs  after  the  application  of  the  Hilbert  transform.  For  all  four 
graphical  comparisons,  the  similarities  and  differences  were  discussed,  as  well  as 
possible  implications  for  application  of  signal  processing. 


118 


V.  Conclusions  and  Future  Work 


The  activities  of  this  research  and  the  results  obtained  were  detailed  as  follows. 
The  comparisons  and  conclusions  drawn  from  the  graphs  generated  were  highlighted. 
Also,  overall  concluding  thoughts  related  to  the  objectives  of  this  evaluation  were 
detailed  and  discussed  to  determine  the  implications  of  the  results.  The  research 
conclusions  were  also  the  basis  for  future  work  concerning  RCS  data  analysis  and 
possible  helpful  future  signal  analysis  tools,  including  the  HHT. 

5.1.  Overall  Summary 

Chapter  I  detailed  the  objectives  of  this  research  used  to  assess  the  effectiveness 
of  using  the  CEMD  method  to  analyze  RCS  data,  as  well  as  compared  the  application  of 
the  Hilbert  transform  and  the  FFT  on  the  decomposed  data  sets  after  the  CEMD 
technique.  To  meet  the  research  objectives,  the  CEMD  algorithm  code  was  found  and 
modified.  The  code  was  also  used  to  apply  the  FFT,  Hilbert  transform,  and  windowed 
FFT  to  produce  graphical  results  for  comparison.  Comparative  analysis  of  the  four 
graphically  generated  results  was  performed.  In  addition,  the  generated  plots  were 
analyzed  to  determine  the  algorithm’s  ability  to  meet  the  objectives  of  this  research 
effort. 

5.2.  Key  Results 

Results  and  conclusions  from  the  work  performed  in  this  thesis  are  presented  as 
follows. 
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1.  Results  for  the  ability  of  the  CEMD  algorithm  to  decompose  the  RCS  real- 
world  data  set. 

Overall,  the  comparison  discussed  in  this  section  displayed  the  functionality  of  the 
CEMD  method  for  the  decomposition  of  a  real-world  data  set.  The  CEMD  method  is  a 
possible  tool  for  decomposition  of  RCS  signals,  when  compared  in  the  time  domain. 

2.  Results  from  the  comparison  of  the  CEMD  algorithm  before  and  after  the 
application  of  the  Elilbert  transform. 

From  the  original  signal,  few  differences  were  noticed  between  the  pre  and  post 
Hilbert  transform  application  magnitude  plots.  No  new  information  could  be  determined 
by  comparison  of  the  plots  represented  by  IMF  #1. 

After  comparing  the  two  IMF  #2  plots,  a  slight  difference  in  the  intensity  of  each 
return  value  was  noticed  with  respect  to  time.  On  the  graph  of  the  IMF  #2  of  the  signal 
before  the  application  of  the  Hilbert  transform,  there  was  approximately  double  the 
intensity  in  time  as  compared  to  the  plot  after  the  application  of  the  Hilbert  transform. 
There  was  a  slight  magnitude  difference,  with  the  primary  and  most  noticeable  difference 
exhibited  by  the  intensity  at  each  time  value  before  the  application  of  the  Hilbert 
transform. 

Beginning  with  IMF  #5,  the  curve  after  the  application  of  the  Hilbert  transform 
showed  a  more  fluid  representation  than  previous  IMFs,  as  well  as  displayed  a  more 
specific  value  representation  of  the  real-world  signal.  After  comparing  the  plots  of  IMF 
#6,  this  plot  showed  that,  after  the  application  of  the  Hilbert  transform,  the  plot  was  more 
fluid  than  before  the  application  of  the  Hilbert  transform.  The  same  characteristics  were 
present  in  both  plots  and,  through  application  of  the  Hilbert  transform,  the  plot  of  IMF  #6 
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was  more  focused  and  displayed  a  cleaner  representation  of  the  real-world  signal  than 
previous  IMF  plots. 

The  signal  decomposition  displayed  in  IMFs  #7  through  #11  showed  the 
sinusoidal  nature  of  the  CEMD  method;  nevertheless,  the  sinusoid  was  more 
distinguishable  after  application  of  the  Hilbert  transform.  Prior  to  the  application  of  the 
Hilbert  transform,  the  decomposition  of  the  signal  demonstrated  noisier  and  contained 
possible  extraneous  signal  values  of  the  decomposed  signal  interpretation. 

3.  Results  from  the  comparison  of  the  FFT  plots  of  the  CEMD  algorithm 
before  and  after  the  application  of  the  Hilbert  transform. 

The  main  difference  between  the  FFT  plot  before  and  after  the  Hilbert  transform 
was  noted  by  the  disappearance  of  the  negative  frequencies  after  the  application  of  the 
Hilbert  transform.  Yet,  before  the  application  of  the  Hilbert  transform,  the  full  frequency 
spectrum  displayed  the  radar  return  signal.  It  was  found  that  the  characteristics  of  the 
graphical  results  occurred  due  to  an  inherent  property  of  the  Hilbert  transform. 

Another  difference  occurred  with  the  magnitude  and  was  indicative  of  the  HHT 
not  being  a  useful  analysis  tool  in  comparison  with  the  FFT  for  the  real-world  data  set 
provided  for  analysis.  A  deeper  understanding  of  exactly  what  the  Hilbert  transform 
application  yields  from  an  RCS  data  set  is  necessary  to  determine  the  usability  of  the 
HHT  in  continued  RCS  data  or  to  determine  whether  the  use  of  the  HHT  is  or  is  not  a 
valid  tool  for  the  type  of  data  used  in  the  intelligence  community. 
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4. 


Results  from  the  comparison  of  the  DTI  plots  of  the  CEMD  algorithm 


before  and  after  the  application  of  the  Hilbert  transform. 

Through  the  comparison  of  the  DTI  plots  both  before  and  after  the  application  of 
the  Hilbert  transform  to  the  original  signal  and  decomposed  data,  it  was  concluded  the 
application  of  the  Hilbert  transform  did  not  provide  new  or  significant  analysis  of  the 
data.  After  the  detailed  comparison  of  each  IMF,  the  Hilbert  transform  application 
seemed  to  delete  approximately  half  of  the  original  frequency  data,  as  well  as  incorrectly 
represented  the  signal  along  the  time  axis. 

5.3.  Concluding  Thoughts 

Various  generated  graphical  results  were  analyzed  and  evaluated  throughout  the 
course  of  this  research  effort.  The  two  principle  objectives  were  investigated.  These  two 
objectives  were:  (1)  to  assess  the  effectiveness  of  the  CEMD  method  and  HHT  as  signal 
analysis  tools  for  complex  RCS  data  and  (2)  to  determine  whether  the  HHT  provided  an 
enhanced  fidelity  and  improved  analysis  of  complex  RCS  data  than  through  use  of  the 
FFT  as  the  analysis  tool.  Through  the  four  comparative  analyses,  these  objectives  were 
met.  The  four  comparisons  performed  were:  (1)  the  comparison  of  the  decomposed  data 
using  the  CEMD  with  various  canonical  shapes,  (2)  the  comparison  of  the  decomposed 
data  using  the  CEMD  with  the  decomposed  data  after  the  application  of  the  Hilbert 
transform,  (3)  the  comparison  of  the  FFT  of  the  decomposed  data  using  the  CEMD  both 
before  and  after  the  Hilbert  transform  was  applied,  and  (4)  the  comparison  of  the  DTI 
plots  of  the  decomposed  data  using  the  CEMD  both  before  and  after  the  application  of 
the  Hilbert  transform. 
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The  first  graphical  comparison  displayed  that  the  CEMD  method  inherent  to  the 
HHT  effectively  decomposed  the  RCS  data  and  exhibited  the  IMFs  of  the  real-world 
signal  in  a  correct  manner.  The  second,  third,  and  fourth  comparisons  were  used  to 
determine  whether  the  hypothesis  that  the  HHT  might  provide  an  enhanced  fidelity  and 
new  or  useful  information  of  the  input  RCS  signal  was  proven.  From  the  comparative 
analyses,  the  application  of  the  Hilbert  transform  did  not  provide  an  enhanced  fidelity 
than  the  application  of  the  FFT.  The  Hilbert  transform  even  appeared  to  remove 
important  data  and  incorrectly  displayed  the  desired  amount  of  fidelity  the  application  of 
the  FFT  provided.  Nevertheless,  it  could  not  be  determined  whether  the  complete  HHT 
was  a  useful  signal  analysis  tool  for  RCS  data;  it  was  concluded  that  the  CEMD  method 
inherent  to  the  HHT  provided  a  useful  analysis  tool  for  RCS  data.  Finally,  compared 
with  the  application  of  the  FFT,  the  Hilbert  transform  inadequately  represented  the  data 
and  was  more  of  a  hindrance  to  the  analysis  of  the  various  plots  generated  of  the  real- 
world  data  set. 

5.4.  Future  Work 

A  list  of  possible  areas  for  future  research  concerning  the  usability  of  the  CEMD 
and  EMD  algorithms,  as  well  as  the  HHT  and  other  mathematical  transforms  as  possible 
analysis  tools,  are  described  as  follows. 

•  Investigation  of  blind  studies  of  the  CEMD  algorithm  should  be  performed.  It  is 
essential  to  establish  the  validity  of  the  CEMD  for  determining  objects  and  the 
shape  of  those  objects  without  prior  knowledge. 
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•  Investigations  of  comparing  various  rocket  bodies  should  be  conducted  to 
distinguish  between  similar  shapes  through  decomposing  and  analyzing  the  IMFs. 
It  is  important  to  know  the  minute  details  for  expedited  analysis  of  the  rocket 
bodies. 

•  Implementation  of  real-valued  RCS  data  as  opposed  to  complex-valued  data 
should  be  used  to  determine  if  similar  results  occur  as  with  the  CEMD  algorithm. 

•  Investigation  of  the  using  time-frequency  transforms  as  possible  new  signal 
analysis  tools  on  decomposed  signals  should  be  examined  compared  to  current 
results  of  use  of  the  HHT  and  FFT. 
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Appendix  A 


Appendix  A  displays  Dr.  Flandrin’s  original  BEMD  algorithm  code,  as  applied  to 
an  example  complex  data  set. 


%bivariate  EMD  illustration ,m 

%illustration  of  the  bivariate  EMD  extension  on  a  real-world 
oceanographic  signal 

%reproduces  Fig.  3  in  "Bivariate  Empirical  Mode  Decomposition",  G. 
Rilling, 

%P.  Flandrin,  P.  Goncalves  and  J.  M.  Lilly,  IEEE  Signal  Processing 
Letters 

Q, 

O 

%G.  Rilling  3/2007  email:  gabriel.rilling@ens-lyon.fr 
load ( ' float_position_record.mat ' , '  x  '  )  ; 

[ imf , nb]  =  cemdc2_f ix ( [ ] , x, 10 ,  [ ] , 32 )  ; 
x  =  hilbert (x) 
n  =  size (imf, 1) ; 

figtitlel  =  'Float  position  record'; 
figure ( ' name ' , figtitlel) 
plot (x) ; 

xlabel (' Displacement  East  (km)  -  Real  part') 

ylabel (' Displacement  North  (km)  -  Imaginary  part') 

title (figtitlel) 
axis  equal; 

set(gca,  'Ylim',  [-250,300]) 

figtitle2  =  'Bivariate  Empirical  Mode  Decomposition  of  Float  signal'; 

figure ( ' name ' , figtitle2) 

subplot  (n+1 , 1 , 1 ) 

plot (real (x) ) 

hold  on 

plot (imag (x) , ' k-- ' ) 

axis  tight 

ylabel ( ' signal ' ) 

title (figtitle2) 

set (gca, ' XTickLabel ' , { } ) 

minmin  =  @ (x) min (x ( : ) ) ; 

maxmax  =  @ (x) max (x ( : ) ) ; 

m  =  minmin ( [real (imf (1 : end-1, :)); imag (imf (1 : end-1, :))]) ; 

M  =  maxmax ( [real ( imf ( 1 : end-1 ,:)); imag ( imf ( 1 : end-1 ,:))]) ; 
for  k  =  1 : n 

subplot (n+1 , 1 , k+1 ) 
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plot (real (imf (k, : ) ) ) 
hold  on 

plot (imag (imf (k, : ) ) , ' k-- ' ) 
axis ( [ 1 , length (x) , m, M] ) 
ylabel ( [ ' d_' , int2str (k) ] ) 
if  k<n 

set (gca, ' XTickLabel ' , { } ) 
end 
end 

ylabel ( ' res .  '  ) 
xlabel('Time  (days) ') 
axis  tight 

Figure  91:  BEMD  Algorithm  by  Dr.  Flandrin 
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Appendix  B 


Appendix  B  shows  the  detailed  code  used  for  producing  the  graphical  results. 


%%  Section  B.l 

%%  Load  all  the  data  and  remake  it  to  work  for  the  algorithm 

load ( ' data .mat ' ) ; 
dataR  =  nb  data.nb  rcs_r'; 
datal  =  nb  data.nb  res  i'; 
prf=round(nb  data.nb  avg_prf ) ; 

Ts  =  round(nb  data.nb  avg_prf ) ; 
data  =  (dataR+i*dataI )  . '  ; 
lengthData=length (data)  ; 

Figure  92:  Section  B.l  Modified  BEMD  Algorithm 


%%  Section  B.2 

%%  Apply  the  CEMD  to  the  inputted  signal 

[imf , nb]  =  cemdc2 ( [ ] , data, [ ] , [ ] ) ; 
n  =  size (imf, 1) ; 

Figure  93:  Section  B.2  Modified  BEMD  Algorithm 


%%  Section  B.3 

%%  Plots  the  Complex  Empirical  Mode  Decomposition  of  the  signal 

figure  ( ) ; 
subplot (2,1,1) 
plot (abs (data) ) 

title ('CEMD  Magnitude  of  Real-World  Signal--Original  Signal'); 

xlabel ( ' Time ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot  (2,1,2) 

plot (unwrap (angle (data) ) ) 

title ('CEMD  Phase  of  Real-World  Signal--Original  Signal'); 

xlabel ( ' Time ' ) 

ylabel (' Angle  (radians) ') 

axis  tight; 

for  k  =  1 : n-1 
figure  ( ) ; 
subplot (2,1,1) 
plot (abs (imf (k, : ) ) ) 

title (['CEMD  Magnitude  of  Real-World  Signal--IMF  #  ' , int2str (k) ] ) ; 

xlabel ( ' Time ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 
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plot (unwrap (angle (imf (k, : ) ) ) ) 

title ( [ ' CEMD  Phase  of  Real-World  Signal--IMF  #  ' , int2str (k) ] ) ; 
xlabel ( ' Time ' ) 
ylabel (' Angle  (radians)  ') 
axis  tight; 
end 

for  k=n : n 
figure  ( ) ; 
subplot (2,1,1) 
plot (abs (imf (k, : ) ) ) 

title ('CEMD  Magnitude  of  Real-World  Signal--Residue ' ) ; 

xlabel ( ' Time ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot  (2,1,2) 

plot (unwrap (angle (imf (k, : ) ) ) ) 

title ('CEMD  Phase  of  Real-World  Signal--Residue ' ) ; 

xlabel ( ' Time ' ) 

ylabel (' Angle  (radians) ') 

axis  tight; 

end 

Figure  94:  Section  B.3  Modified  BEMD  Algorithm 


%%  Section  B.4 

%%  Take  the  HT  of  the  Decomposed  data  and  plot  the  HT  of  the  CEMD  of 
the 

%%  data 

figure  ( ) ; 

hilbertOriginal=hilbert (data)  ; 
hilbertMagOriginal=abs (hilbert (real (data) ) ) ; 
hilbertPhaseOriginal=unwrap (angle (hilbert (imag (data) ) ) )  ; 
subplot (2,1,1) 
plot  (hilbertMagOriginal ) 

title (' Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal--Original 

Signal ' ) ; 

xlabel ( ' Time ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot  (2,1,2) 

plot (hilbertPhaseOriginal ) 

title (' Hilbert-Huang  Transform  Phase  of  Real-World  Signal--Original 

Signal ' )  ; 

xlabel ( ' Time ' ) 

ylabel (' Angle  (radians) ') 

axis  tight; 

for  k  =  1 : n-1 

hilbertMag (k,  : ) =abs (hilbert (real (imf (k,  : ) ) ) )  ; 
hilbertPhase (k,  : ) =unwrap (angle (hilbert (imag (imf (k,  :))))); 
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figure  () 
subplot (2,1,1) 
plot (hilbertMag (k, : ) ) 

title ([' Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal--IMF 
#  ' ,  int2str (k) ] ) 

ylabel ( 'RCS  (dBsm) ' ) 
xlabel (' Time ') axis  tight; 
subplot (2,1,2) 
plot  (hilbertPhase (k,  : ) ) 

title ([' Hilbert-Huang  Transform  Phase  of  Real-World  Signal--IMF  #  ', 
int2str (k) ] ) 

ylabel (' Angle  (radians) ') 
xlabel ( ' Time ' ) 
axis  tight; 

end 

for  k  =  n:n 

hilbertMag (k, : ) =abs (hilbert (real (imf (k, : ) ) ) ) ; 

hilbertPhase (k, : ) =unwrap (angle (hilbert (imag (imf (k, :))))); 

figure  () 

subplot (2,1,1) 

plot (hilbertMag (k, : ) ) 

title (' Hilbert-Huang  Transform  Magnitude  of  Real-World  Signal-- 
Residue ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

xlabel ( ' Time ' ) 

axis  tight; 

subplot  (2 , 1,2) 

plot (hilbertPhase (k, : ) ) 

title  (' Hilbert-Huang  Transform  Phase  of  Real-World  Signal--Residue ' ) 
ylabel (' Angle  (radians) ') 
xlabel ( ' Time ' ) 
axis  tight; 

end 

Figure  95:  Section  B.4  Modified  BEMD  Algorithm 


%%  Section  B.5 

%%  Take  the  FFT  of  the  decomposed  signal  before  HT  and  plot  the  graphs 
figure  ( ) ; 

fftOriginal= (1/n) . *fft (data) ; 
f f tMagOriginal=abs ( f f tOriginal ) ; 
f f tPhaseOriginal=unwrap (angle ( ff tOriginal ) ) ; 
subplot (2,1,1) 

plot (fftshift ( f f tMagOriginal ) ) 

title ('FFT  Magnitude  of  Real-World  Signal--Original  Signal'); 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot ( f f tPhaseOriginal ) 

title ('FFT  Phase  of  Real-World  Signal--Original  Signal'); 
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xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians) ') 
axis  tight; 

for  k  =  1 : n-1 
figure ( ) ; 

f f tCEMD ( k, : ) = (1/n) . *f ft (imf (k, :)); 
fftMagCEMD (k,  : ) =abs (fftCEMD (k,  : ) )  ; 
fftPhaseCEMD (k,  : ) =unwrap (angle (imf (k,  : ) ) ) ; 
subplot (2,1,1) 

plot (fftshift (fftMagCEMD (k,  : ) )  ) 

title (('FFT  Magnitude  of  Real-World  Signal--IMF  #  '  , int2str (k) ] ) ; 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot (fftPhaseCEMD (k,  :  )  ) 

title (['FFT  Phase  of  Real-World  Signal--IMF  #  '  , int2str (k) ] )  ; 
xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians) ') 
axis  tight; 
end 

for  k=n : n 
figure  ( ) ; 

fftCEMD (k,  :  )  =  (1/n) *fft (imf (k,  :)); 
fftMagCEMD (k, : ) =abs (fftCEMD (k, : ) ) ; 
fftPhaseCEMD (k,  : ) =unwrap (angle (imf (k,  : ) ) ) ; 
subplot (2,1,1) 

plot  (fftshift (fftMagCEMD (k,  :  )  )  ) 

title (' FFT  Magnitude  of  Real-World  Signal--Residue ' ) ; 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot (fftPhaseCEMD (k, : ) ) 

title ('FFT  Phase  of  Real-World  Signal--Residue ' ) ; 
xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians) ') 
axis  tight; 
end 


Figure  96:  Section  B.5  Modified  BEMD  Algorithm 


%%  Section  B.6 

%%  Take  the  FFT  of  the  decomposed  signal  after  HT  and  plot  the  graphs 
figure  ( ) ; 

f f tOriginalHilbert=f f t (hilbert (data)  )  ; 

f f tMagOriginalHilbert=abs ( (1/n)  . *fft (hilbert (real (data) ) ) )  ; 
f f tPhaseOriginalHilbert=unwrap (angle ( (hilbert (imag (data) ) ) ) ) ; 
subplot (2,1,1) 
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plot ( (fftshift (fftMagOriginalHilbert) ) ) 

title ( ' HHT  of  FFT  Magnitude  of  Real-World  Signal--Original  Signal'); 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot (fftPhaseOriginalHilbert) 

title ('HHT  of  FFT  Phase  of  Real-World  Signal--Original  Signal'); 
xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians) ') 
axis  tight; 

for  k  =  1 : n-1 
figure ( ) ; 

f f tMagHilbert (k, : ) =abs ( (1/n) . *fft (hilbert (real (imf (k, :))))) ; 
f f tPhaseHilbert (k,  : ) =unwrap (angle (fft (hilbert (imag (imf (k,  :))))))  ; 
subplot (2,1,1) 

plot (fftshift ( ff tMagHilbert (k, : ) ) ) 

title (('HHT  of  FFT  Magnitude  of  Real-World  Signal--IMF  # 

'  ,  int2str (k)  ]  )  ; 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot ( ff tPhaseHilbert (k,  : )  ) 

title (['HHT  of  FFT  Phase  of  Real-World  Signal--IMF  #  '  , int2str (k) ] )  ; 
xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians)  ') 
axis  tight; 
end 

for  k=n : n 
figure  ( ) ; 

ff tMagHilbert (k,  : ) =abs (  (1/n)  . *fft (hilbert (real (imf  (k,  :))))) ; 
ff tPhaseHilbert (k,  : ) =unwrap (angle (fft (hilbert (imag (imf (k,  :))))))  ; 
subplot  (2,1,1) 

plot (fftshift ( ff tMagHilbert (k, : ) ) ) 

title ('HHT  of  FFT  Magnitude  of  Real-World  Signal--Residue ' )  ; 

xlabel ( ' Frequency ' ) 

ylabel ( 'RCS  (dBsm) ' ) 

axis  tight; 

subplot (2,1,2) 

plot ( ff tPhaseHilbert (k,  : )  ) 

title  ('HHT  of  FFT  Phase  of  Real-World  Signal--Residue ' ) ; 
xlabel ( ' Frequency ' ) 
ylabel (' Angle  (radians)') 
axis  tight; 
end 

Figure  97:  Section  B.6  Modified  BEMD  Algorithm 
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%%  Section  B.7 

%%  Take  the  windowed  FFT  of  the  decomposed  signal  before  HT  to  generate 
the 

%%  DTI  plots 

%Define  resolution  for  the  windowed  FFT 

step_res=50 ; 

sec=96; 

%Calculate  and  store  the  windowed  FFT  value  of  the  original  signal 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

original (b, : ) =f ft ( [data (counter*step_res+l : counter*step_res+l+prf ) , . . . 
zeros (1, 1024-prf ) ] ) ; 
counter=counter+l  ; 
end 

%Graph  the  windowed  FFT  of  the  original  signal 
y=linspace ( -50 , 50 ,  1025); 

x=linspace ( 0 ,  100,  prf* (sec-1) /step_res+l) ; 
figure  ( ) 

pcolor (x,  y,  fftshift (20*log (abs (original) ') ,  1)); 
shading  interp 

title ('DTI  Plot--CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-- 

Original  Signal') 

xlabel ( ' Time ' ) 

ylabel (' Frequency  (GHz) ') 

axis  tight; 

caxis ( [-40  40] )  ; 

colorbar ; 

%Get  the  windowed  FFT  of  each  of  the  IMFs  and  plot  the  windowed  FFT  of 
each 
%  IMF 

for  a  =  1 : n-1 

sig=imf (a, : ) ; 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

ff(b, :)  =  f ft ( [sig ( counter* step_res+l : counter *step_res+l+prf) , . . . 

zeros (1, 1024-prf) ] ) ; 
counter=counter+l  ; 
end 

figure ( ) 

pcolor (x,y,  fftshift(20*log(abs(ff)  ') ,1) ) ; 
shading  interp; 

title (['DTI  Plot--CEMD  Magnitude  of  the  FFT  of  Real-World  Signal--IMF 

#  '  .  .  . 

,  int2str (a) ] ) 
xlabel ( ' Time ' ) 
ylabel (' Frequency  (GHz)') 
axis  tight; 
caxis ( [-40  40] ); 
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colorbar; 

end 


%Windowed  FFT  of  the  residue  of  the  signal 
for  a  =  n : n 

sig=imf (a, : ) ; 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

ff(b, :)  =  f ft ( [sig ( counter* step_res+l : counter *step_res+l+prf) , . . . 

zeros ( 1 , 1024-prf ) ] ) ; 
counter=counter+l  ; 
end 

figure  ( ) 

pcolor(x,y,  fftshift(20*log(abs(ff)  ') ,1) ) ; 
shading  interp; 

title ('DTI  Plot--CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-- 
Residue  ' ) 
xlabel ( ' Time ' ) 
ylabel (' Frequency  (GHz)') 
axis  tight; 
caxis (  [-4  0  40]); 
colorbar; 
end 

Figure  98:  Section  B.7  Modified  BEMD  Algorithm 


%%  Section  B.8 

%%  Take  the  windowed  FFT  of  the  decomposed  signal  after  HT  to  generate 
the  DTI  plots 

%Define  resolution  for  the  fft 

step_res=50 ; 

sec=96; 

%Calculate  and  store  the  windowed  FFT  value  of  the  original  signal 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

originalH (b, : ) =fft ( [hilbert (data (counter*step_res+l : . . . 

counter*step_res+l+prf ) ) , zeros (1, 1024-prf) ] ) ; 
counter=counter+l ; 

end 

%Graph  the  windowed  FFT  of  the  original  signal 
fftMagOriginalHilbert=f ft shift (20*log(abs (originalH' ) ) ) ; 
y=linspace ( -50 , 50 ,  1025); 

x=linspace ( 0 ,  100,  prf* (sec-1) /step_res  +  l) ; 
figure  ( ) 

pcolor (x, y, f ftMagOriginalHilbert) 
shading  interp 

title ('DTI  Plot--HHT  Magnitude  of  the  FFT  Real-World  Signal--Original 
Signal ' ) 
xlabel ( ' Time ' ) 
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ylabel (' Frequency  (GHz)') 
axis  tight; 
caxis ( [-40  40] )  ; 
colorbar; 

%Get  the  windowed  FFT  of  each  of  the  IMFs  after  the  HT  has  been  applied 
to  the  signal  and  plot  the  windowed  FFT  of  each  IMF 
for  a  =  1 : n-1 

sigH=hilbert (imf (a, : ) ) ; 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

ffH(b, :)  =  fft ( [sigH (counter*step  res+1 : counter*step  res+l+prf ) , . . . 

zeros (1, 1024-prf ) ] ) ; 
counter=counter+l ; 

end 

figure ( ) 

fftDTIMagHilbert=f ft shift (20*log (abs (ffH ' ) ) ) ; 
pcolor(x,y,  fftDTIMagHilbert); 
shading  interp; 

title (['DTI  Plot--HHT  Magnitude  of  the  FFT  Real-World  Signal--IMF  # 

I 

,  int2str (a) ] ) 
xlabel ( ' Time ' ) 
ylabel (' Frequency  (GHz)') 
axis  tight; 
caxis ( [-40  40] )  ; 
colorbar; 
end 

%Windowed  FFT  of  the  residue  of  the  signal 
for  a  =  n:n 

sigH=hilbert (imf (a,  : )  )  ; 
counter=0 ; 

for  b=l :prf* (sec-1) /step  res+1 

ffH (b, : )  =  f ft ( [sigH ( counter* step_res+l : counter *step_res+l+prf) , . . . 

zeros ( 1 , 1024-prf) ] ) ; 
counter=counter+l ; 
end 

figure ( ) 

fftDTIMagHilbert=f ft shift (20*log (abs (ffH ' ) ) ) ; 
pcolor(x,y,  fftDTIMagHilbert); 
shading  interp; 

title ('DTI  Plot--HHT  Magnitude  of  the  FFT  Real-World  Signal--Residue ' ) 

xlabel ( ' Time ' ) 

ylabel (' Frequency  (GHz)') 

axis  tight; 

caxis ( [-40  40] )  ; 

colorbar; 

end 

Figure  99:  Section  B.7  Modified  BEMD  Algorithm 
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Appendix  C 


Appendix  C  explains  the  inputs  and  outputs  of  cemdc  f  ix  ( )  and  explains  what 
task  each  part  of  the  function  performs. 


^CEMDC_FIX  bivariate  Empirical  Mode  Decomposition,  first  algorithm 


Syntax 


%  [ IMF, NB  ITERATIONS] =CEMDC  FIX(T,X,NB  ITERATIONS , MAX  IMFS,NDIRS); 


Description 


%  computes  bivariate  EMD,  first  algorithm  [1]  with  NB^ITERATONS  sifting 
%  iterations  for  each  IMF 

O, 

o 

%  mean  of  boolean  array  { (mean  amplitude) / (envelope_amplitude)  > 
THRESHOLD}  <  TOLERANCE 


%  inputs: 


sampled . 


to 


T:  sampling  times.  If  T= [ ]  ,  the  signal  is  assumed  uniformly 
X:  analyzed  signal 

NB  ITERATIONS:  number  of  sifting  iterations  to  be  performed 


%  extract  each  IMF.  If  NB  ITERATIONS  is  empty  or  unspecified, 

10  iterations 

%  are  performed  by  default. 

%  Note:  The  effective  number  of  sifting  iterations  might  be 


less 

g, 
o 

has 

g, 
o 

g, 
o 

is 

g. 
o 

extract  as 


than  NB  ITERATIONS  for  the  last  modes  if  the  sifting  process 
to  be  stopped  because  of  a  lack  of  extrema. 

MAX  IMFS :  maximum  number  of  IMFs  to  be  extracted.  If  MAX  IMFS 
zero,  empty  or  unspecified,  the  default  behavior  is  to 
many  IMFs  as  possible. 

NDIRS:  number  of  directions  used  to  compute  the  local  mean. 

If  unspecified,  the  default  value  is  4. 

rem:  the  actual  number  of  directions  (according  to  [1])  is 


2  *NDIRS 

Q. 

O 

%  outputs: 

g,  _ 

o 

o, 

o 

each  mode 


IMF:  intrinsic  mode  functions  (IMFs)  (last  line  =  residual) 
NB  ITERATIONS:  effective  number  of  sifting  iterations  for 
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Examples : 


g, 

o 

o, 

o 

o, 

o 

o, 

o 


%  workspace  %  T:  lxN  time  instants 
%  X:  lxN  signal  data 


o, 

o 


%»IMF  =  CEMDC_  

%»  [  IMF,  NB_IT]  =  CEMDC_FIX([]  

%»IMF  =  CEMDC_  

%»  [  IMF,  NB_IT ]  =  CEMDC  FIX  (  [  ]  ,  X,  [  ]  ,  4  )  ; 


O, 

o 

o, 

o 


References 


o, 

o 

o, 

o 


%  [1]  G.  Rilling,  P.  Flandrin,  P.  Gongalves  and  J.  M.  Lilly., 
%  "Bivariate  Empirical  Mode  Decomposition", 

%  Signal  Processing  Letters  (submitted) 


O, 


o 


o, 

o 


%  See  also 

%  (c)emd  visu  (visualization), 

%  emd  (slow  but  has  many  options), 

%  cemdc,  cemdc2,  cemdc2_fix  (other  fast  implementations  of  bivariate 
EMD) 


O, 


o 


o. 


o 


%  G.  Rilling,  last  modification:  3.2007 
%  gabriel.rilling@ens-lyon.fr 


g. 

o 


%  code  based  on  a  student  project  by  T.  Boustane  and  G.  Quellec, 
11.03.2004 

%  supervised  by  P.  Chainais  (ISIMA  -  LIMOS  -  Universite  Blaise  Pascal  - 
Clermont  II 

%  email  :  pchainai@isima.fr) . 


Figure  100:  CEMD_FIX  function 
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Appendix  D 


The 

application. 


IMF  #3  and  #4  plots  of  DTI  before  and  after  the  Hilbert  transform 


DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #3 
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Figure  101:  DTI  Magnitude  of  IMF  #3 

DTI  Plot-HHT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #3 
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Figure  102:  DTI  HHT  Magnitude  of  IMF  #3 
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DTI  Plot-CEMD  Magnitude  of  the  FFT  of  Real-World  Signal-IMF  #4 


Time  (s) 


Figure  103:  DTI  Magnitude  of  IMF  #4 

DTI  Plot- FI  FIT  Magnitude  of  the  FFT  Real-World  Signal-IMF  #4 


Time  (s) 

Figure  104:  DTI  HHT  Magnitude  of  IMF  #4 
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